From 145bf00ef77abfd838cb0de25b54b9f88fe18b50 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 30 Aug 2017 22:42:39 -0700 Subject: [PATCH 01/34] use properties for n, x0 for the base mesh --- discretize/BaseMesh.py | 84 +++++++++++++++++++++++------------------- 1 file changed, 47 insertions(+), 37 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index b35e309ef..468fee005 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -1,42 +1,52 @@ import numpy as np +import properties from discretize import utils -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): - - # Check inputs - if x0 is None: - x0 = np.zeros(len(n)) - - if not len(n) == len(x0): + # Properties + n = properties.Array( + "number of cells in each direction (dim, )", + dtype=int, + required=True, + shape=('*',) + ) + + x0 = properties.Array( + "origin of the mesh (dim, )", + dtype=float, + shape=('*',), + required=True + ) + + def __init__(self, n, **kwargs): + self.n = n # number of dimensions + super(BaseMesh, self).__init__(**kwargs) + + # set x0 to zeros if it is none + if self.x0 is None: + self.x0 = np.zeros(len(self.n)) + + # Validators + @properties.validator('n') + def check_n_shape(self, change): + 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']) + ) + + @properties.validator('x0') + def check_x0_vs_n(self, change): + if len(self.n) != len(change['value']): raise Exception("Dimension mismatch. x0 != len(n)") - if len(n) > 3: - raise Exception("Dimensions higher than 3 are not supported.") - - # 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) - - @property - def x0(self): - """Origin of the mesh - - :rtype: numpy.array - :return: x0, (dim, ) - """ - return self._x0 - @property def dim(self): """The dimension of the mesh (1, 2, or 3). @@ -44,7 +54,7 @@ def dim(self): :rtype: int :return: dim """ - return self._dim + return len(self.n) @property def nC(self): @@ -293,8 +303,8 @@ def projectEdgeVector(self, eV): class BaseRectangularMesh(BaseMesh): """BaseRectangularMesh""" - def __init__(self, n, x0=None): - BaseMesh.__init__(self, n, x0) + def __init__(self, n, **kwargs): + BaseMesh.__init__(self, n, **kwargs) @property def nCx(self): @@ -303,7 +313,7 @@ def nCx(self): :rtype: int :return: nCx """ - return int(self._n[0]) + return int(self.n[0]) @property def nCy(self): @@ -314,7 +324,7 @@ def nCy(self): """ if self.dim < 2: return None - return int(self._n[1]) + return int(self.n[1]) @property def nCz(self): @@ -325,7 +335,7 @@ def nCz(self): """ if self.dim < 3: return None - return int(self._n[2]) + return int(self.n[2]) @property def vnC(self): @@ -622,7 +632,7 @@ def outKernal(xx, nn): def switchKernal(xx): """Switches over the different options.""" if xType in ['CC', 'N']: - nn = (self._n) if xType == 'CC' else (self._n+1) + nn = (self.n) if xType == 'CC' else (self._n+1) assert xx.size == np.prod(nn), "Number of elements must not change." return outKernal(xx, nn) elif xType in ['F', 'E']: From 19882a238977240b615a73735fcc1bbb7898803b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 30 Aug 2017 22:44:03 -0700 Subject: [PATCH 02/34] use n instead of private variable --- discretize/BaseMesh.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index 468fee005..5a5c012c3 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -71,7 +71,7 @@ def nC(self): M = discretize.TensorMesh([np.ones(n) for n in [2,3]]) M.plotGrid(centers=True, showIt=True) """ - return int(self._n.prod()) + return int(self.n.prod()) @property def nN(self): @@ -88,7 +88,7 @@ def nN(self): M = discretize.TensorMesh([np.ones(n) for n in [2,3]]) M.plotGrid(nodes=True, showIt=True) """ - return int((self._n+1).prod()) + return int((self.n+1).prod()) @property def nEx(self): @@ -97,7 +97,7 @@ def nEx(self): :rtype: int :return: nEx """ - return int((self._n + np.r_[0, 1, 1][:self.dim]).prod()) + return int((self.n + np.r_[0, 1, 1][:self.dim]).prod()) @property def nEy(self): @@ -108,7 +108,7 @@ def nEy(self): """ if self.dim < 2: return None - return int((self._n + np.r_[1, 0, 1][:self.dim]).prod()) + return int((self.n + np.r_[1, 0, 1][:self.dim]).prod()) @property def nEz(self): @@ -119,7 +119,7 @@ def nEz(self): """ if self.dim < 3: return None - return int((self._n + np.r_[1, 1, 0][:self.dim]).prod()) + return int((self.n + np.r_[1, 1, 0][:self.dim]).prod()) @property def vnE(self): @@ -158,7 +158,7 @@ def nFx(self): :rtype: int :return: nFx """ - return int((self._n + np.r_[1, 0, 0][:self.dim]).prod()) + return int((self.n + np.r_[1, 0, 0][:self.dim]).prod()) @property def nFy(self): @@ -169,7 +169,7 @@ def nFy(self): """ if self.dim < 2: return None - return int((self._n + np.r_[0, 1, 0][:self.dim]).prod()) + return int((self.n + np.r_[0, 1, 0][:self.dim]).prod()) @property def nFz(self): @@ -180,7 +180,7 @@ def nFz(self): """ if self.dim < 3: return None - return int((self._n + np.r_[0, 0, 1][:self.dim]).prod()) + return int((self.n + np.r_[0, 0, 1][:self.dim]).prod()) @property def vnF(self): @@ -632,7 +632,7 @@ def outKernal(xx, nn): def switchKernal(xx): """Switches over the different options.""" if xType in ['CC', 'N']: - nn = (self.n) if xType == 'CC' else (self._n+1) + nn = (self.n) if xType == 'CC' else (self.n+1) assert xx.size == np.prod(nn), "Number of elements must not change." return outKernal(xx, nn) elif xType in ['F', 'E']: From af1035d3b3e9cd3d0d2eef21757f58281a56c683 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 30 Aug 2017 22:50:47 -0700 Subject: [PATCH 03/34] whitespace cleanup in base mesh --- discretize/BaseMesh.py | 63 ++++++++++++++++++++++++++++++++---------- 1 file changed, 49 insertions(+), 14 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index 5a5c012c3..7fca3fca9 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -24,6 +24,7 @@ class BaseMesh(properties.HasProperties): required=True ) + # Instantiate the class def __init__(self, n, **kwargs): self.n = n # number of dimensions super(BaseMesh, self).__init__(**kwargs) @@ -600,16 +601,34 @@ 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 = [ + '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): - 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 @@ -620,7 +639,11 @@ def r(self, x, xType='CC', outType='CC', format='V'): 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.""" @@ -633,10 +656,13 @@ 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' + # 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 nn = self.vnF if xType == 'F' else self.vnE nn = np.r_[0, nn] @@ -648,13 +674,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: @@ -668,7 +701,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]): From 693c5d021eb39da8a89d5108ac9f10eef01eef4f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 30 Aug 2017 22:56:02 -0700 Subject: [PATCH 04/34] explicitly add x0 as an input to the basemesh --- discretize/BaseMesh.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index 7fca3fca9..cfd9faca7 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -25,8 +25,10 @@ class BaseMesh(properties.HasProperties): ) # Instantiate the class - def __init__(self, n, **kwargs): + def __init__(self, n, x0=None, **kwargs): self.n = n # number of dimensions + if x0 is not None: + self.x0 = x0 super(BaseMesh, self).__init__(**kwargs) # set x0 to zeros if it is none @@ -304,8 +306,8 @@ def projectEdgeVector(self, eV): class BaseRectangularMesh(BaseMesh): """BaseRectangularMesh""" - def __init__(self, n, **kwargs): - BaseMesh.__init__(self, n, **kwargs) + def __init__(self, n, x0=None, **kwargs): + BaseMesh.__init__(self, n, x0=None, **kwargs) @property def nCx(self): From 809a149922b470c2003e1bef42e9d2027916ff2b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 30 Aug 2017 23:29:23 -0700 Subject: [PATCH 05/34] esure x0 being passed in init of BaseRectangularMesh --- discretize/BaseMesh.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index cfd9faca7..3e5a2470a 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -25,15 +25,14 @@ class BaseMesh(properties.HasProperties): ) # Instantiate the class - def __init__(self, n, x0=None, **kwargs): + def __init__(self, n, x0=None): self.n = n # number of dimensions - if x0 is not None: - self.x0 = x0 - super(BaseMesh, self).__init__(**kwargs) - # set x0 to zeros if it is none - if self.x0 is None: + # origin of the mesh + if x0 is None: self.x0 = np.zeros(len(self.n)) + else: + self.x0 = x0 # Validators @properties.validator('n') @@ -306,8 +305,8 @@ def projectEdgeVector(self, eV): class BaseRectangularMesh(BaseMesh): """BaseRectangularMesh""" - def __init__(self, n, x0=None, **kwargs): - BaseMesh.__init__(self, n, x0=None, **kwargs) + def __init__(self, n, x0=None): + BaseMesh.__init__(self, n, x0=x0) @property def nCx(self): From 44fea79fd970c1f2ff66fa82118aa9a7ff58c3a6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 30 Aug 2017 23:29:47 -0700 Subject: [PATCH 06/34] minimal injection of properties to base tensor mesh --- discretize/TensorMesh.py | 43 ++++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/discretize/TensorMesh.py b/discretize/TensorMesh.py index 36451da73..4f94c0ae5 100644 --- a/discretize/TensorMesh.py +++ b/discretize/TensorMesh.py @@ -1,6 +1,7 @@ from __future__ import print_function import numpy as np import scipy.sparse as sp +import properties from discretize import utils @@ -18,7 +19,22 @@ class BaseTensorMesh(BaseMesh): _unitDimensions = [1, 1, 1] - def __init__(self, h_in, x0_in=None): + # properties + h = properties.List( + "h is a list containing the cell widths of the tensor mesh in each " + "dimension.", + properties.Array( + "widths of the tensor mesh in a single dimension", + dtype=float, + shape=("*",) + ), + max_length=3 + ) + + def __init__(self, h, x0=None): + h_in = h + x0_in = x0 + assert type(h_in) in [list, tuple], 'h_in must be a list' assert len(h_in) in [1, 2, 3], 'h_in must be of dimension 1, 2, or 3' h = list(range(len(h_in))) @@ -55,38 +71,27 @@ def __init__(self, h_in, x0_in=None): "'C' to center, or 'N' to be negative.".format(i) ) - if isinstance(self, BaseRectangularMesh): - BaseRectangularMesh.__init__( - self, np.array([x.size for x in h]), x0 - ) - else: - BaseMesh.__init__(self, np.array([x.size for x in h]), x0) + super(BaseTensorMesh, self).__init__( + np.array([x.size for x in h]), x0=x0 + ) # Ensure h contains 1D vectors - self._h = [utils.mkvc(x.astype(float)) for x in h] - - @property - def h(self): - """ - h is a list containing the cell widths of the tensor mesh in each - dimension. - """ - return self._h + self.h = [utils.mkvc(x.astype(float)) for x in h] @property def hx(self): "Width of cells in the x direction" - return self._h[0] + return self.h[0] @property def hy(self): "Width of cells in the y direction" - return None if self.dim < 2 else self._h[1] + return None if self.dim < 2 else self.h[1] @property def hz(self): "Width of cells in the z direction" - return None if self.dim < 3 else self._h[2] + return None if self.dim < 3 else self.h[2] @property def vectorNx(self): From d688b8565885543991d9c391acb33c49a8cacc37 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 07:53:19 -0700 Subject: [PATCH 07/34] have tree mesh look at top level h rather than a hidden property --- discretize/TreeMesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index 4e6172a1e..8c7be0397 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -105,8 +105,8 @@ def __init__(self, h, x0=None, levels=None): BaseTensorMesh.__init__(self, h, x0) if levels is None: - levels = int(np.log2(len(self._h[0]))) - assert np.all(len(_) == 2**levels for _ in self._h), "must make h and levels match" + levels = int(np.log2(len(self.h[0]))) + assert np.all(len(_) == 2**levels for _ in self.h), "must make h and levels match" self._levels = levels self._levelBits = int(np.ceil(np.sqrt(levels)))+1 From 209856a57cf1b55f9a35d63429c48488c6baa1d5 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 08:36:34 -0700 Subject: [PATCH 08/34] add ability to save and load tensor meshes --- discretize/BaseMesh.py | 21 ++++++++++++- discretize/TensorMesh.py | 59 ++++++++++++++++++++++------------- discretize/utils/__init__.py | 18 ++++++----- discretize/utils/meshutils.py | 18 +++++++++++ tests/base/test_properties.py | 39 +++++++++++++++++++++++ 5 files changed, 126 insertions(+), 29 deletions(-) create mode 100644 tests/base/test_properties.py diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index 3e5a2470a..d556ee21a 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -1,5 +1,7 @@ import numpy as np import properties +import os +import json from discretize import utils @@ -25,10 +27,11 @@ class BaseMesh(properties.HasProperties): ) # Instantiate the class - def __init__(self, n, x0=None): + def __init__(self, n, **kwargs): self.n = n # number of dimensions # origin of the mesh + x0 = kwargs.pop('x0') if x0 is None: self.x0 = np.zeros(len(self.n)) else: @@ -302,6 +305,22 @@ 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): + """ + 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 + class BaseRectangularMesh(BaseMesh): """BaseRectangularMesh""" diff --git a/discretize/TensorMesh.py b/discretize/TensorMesh.py index 4f94c0ae5..d0f62019e 100644 --- a/discretize/TensorMesh.py +++ b/discretize/TensorMesh.py @@ -31,29 +31,46 @@ class BaseTensorMesh(BaseMesh): max_length=3 ) - def __init__(self, h, x0=None): - h_in = h - x0_in = x0 - - assert type(h_in) in [list, tuple], 'h_in must be a list' - assert len(h_in) in [1, 2, 3], 'h_in must be of dimension 1, 2, or 3' - h = list(range(len(h_in))) - for i, h_i in enumerate(h_in): - if utils.isScalar(h_i) and type(h_i) is not np.ndarray: - # This gives you something over the unit cube. - h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) - elif type(h_i) is list: - h_i = utils.meshTensor(h_i) - assert isinstance(h_i, np.ndarray), ( - "h[{0:d}] is not a numpy array.".format(i) + def __init__(self, h_in=None, **kwargs): + + # Cell widths + if h_in is not None: + # Sanity Checks + assert 'h' not in kwargs.keys(), ( + 'Only 1 descriptor of cell sizes can be used to instantiate ' + 'the class' ) - assert len(h_i.shape) == 1, ( - "h[{0:d}] must be a 1D numpy array.".format(i) + assert type(h_in) in [list, tuple], 'h_in must be a list' + assert len(h_in) in [1, 2, 3], ( + 'h_in must be of dimension 1, 2, or 3' ) - h[i] = h_i[:] # make a copy. + # build h + h = list(range(len(h_in))) + for i, h_i in enumerate(h_in): + if utils.isScalar(h_i) and type(h_i) is not np.ndarray: + # This gives you something over the unit cube. + h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) + elif type(h_i) is list: + h_i = utils.meshTensor(h_i) + assert isinstance(h_i, np.ndarray), ( + "h[{0:d}] is not a numpy array.".format(i) + ) + assert len(h_i.shape) == 1, ( + "h[{0:d}] must be a 1D numpy array.".format(i) + ) + h[i] = h_i[:] # make a copy. + + else: + assert 'h' in kwargs.keys(), 'cell widths must be provided' + h = kwargs.pop('h') + + # Origin of the mesh x0 = np.zeros(len(h)) - if x0_in is not None: + + if 'x0' in kwargs.keys(): + x0_in = kwargs.pop('x0') + assert len(h) == len(x0_in), "Dimension mismatch. x0 != len(h)" for i in range(len(h)): x_i, h_i = x0_in[i], h[i] @@ -526,8 +543,8 @@ class TensorMesh(BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, _meshType = 'TENSOR' - def __init__(self, h_in, x0=None): - BaseTensorMesh.__init__(self, h_in, x0) + def __init__(self, h_in=None, **kwargs): + BaseTensorMesh.__init__(self, h_in, **kwargs) def __str__(self): outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) diff --git a/discretize/utils/__init__.py b/discretize/utils/__init__.py index ed318e7fb..d1c296b30 100644 --- a/discretize/utils/__init__.py +++ b/discretize/utils/__init__.py @@ -1,13 +1,17 @@ from __future__ import print_function -from .matutils import (mkvc, sdiag, sdInv, speye, kron3, spzeros, ddx, av, - av_extrap, ndgrid, ind2sub, sub2ind, getSubArray, - inv3X3BlockDiagonal, inv2X2BlockDiagonal, TensorType, - makePropertyTensor, invPropertyTensor, Zero, - Identity) +from .matutils import ( + mkvc, sdiag, sdInv, speye, kron3, spzeros, ddx, av, + av_extrap, ndgrid, ind2sub, sub2ind, getSubArray, + inv3X3BlockDiagonal, inv2X2BlockDiagonal, TensorType, + makePropertyTensor, invPropertyTensor, Zero, + Identity +) from .codeutils import (isScalar, asArray_N_x_Dim) -from .meshutils import (exampleLrmGrid, meshTensor, closestPoints, - ExtractCoreMesh, random_model) +from .meshutils import ( + load_mesh, exampleLrmGrid, meshTensor, closestPoints, ExtractCoreMesh, + random_model +) from .curvutils import volTetra, faceInfo, indexCube from .interputils import interpmat from .coordutils import rotatePointsFromNormals, rotationMatrixFromNormals diff --git a/discretize/utils/meshutils.py b/discretize/utils/meshutils.py index 5de236c98..c3b6691d0 100644 --- a/discretize/utils/meshutils.py +++ b/discretize/utils/meshutils.py @@ -1,4 +1,6 @@ import numpy as np +import properties +import json from .matutils import ndgrid from .codeutils import asArray_N_x_Dim from .codeutils import isScalar @@ -12,6 +14,22 @@ num_types = [int, float] +def load_mesh(filename): + """ + Open a json file and load the mesh into the target class + + As long as there are no namespace conflicts, the target __class__ + will be stored on the properties.HasProperties registry and may be + fetched from there. + + :param str filename: name of file to read in + """ + with open(filename, 'r') as outfile: + jsondict = json.load(outfile) + data = properties.HasProperties.deserialize(jsondict, trusted=True) + return data + + def exampleLrmGrid(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" diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py new file mode 100644 index 000000000..8d00027cc --- /dev/null +++ b/tests/base/test_properties.py @@ -0,0 +1,39 @@ +import unittest +import os +import discretize + + +class TensorTest(unittest.TestCase): + + n = [4, 5, 9] + x0 = [-0.5, -0.25, 0] + + def setUp(self): + self.mesh = discretize.TensorMesh(self.n, x0=self.x0) + + def test_save_load(self): + print('\nTesting save / load of Tensor Mesh ...') + mesh0 = self.mesh + f = mesh0.save() + mesh1 = discretize.utils.load_mesh(f) + + assert mesh0.nC == mesh1.nC , ( + 'Number of cells not the same, {} != {}'.format(mesh0.nC, mesh1.nC) + ) + + assert (mesh0.x0 == mesh1.x0).all(), ( + 'x0 different. {} != {}'.format(mesh0.x0, mesh1.x0) + ) + + for i in range(len(mesh0.h)): + assert (mesh0.h[i] == mesh1.h[i]).all(), ( + 'mesh h[{}] different'.format(i) + ) + + os.remove(f) + + + + +if __name__ == '__main__': + unittest.main() From 85a4ba00bba0733457dc62e37479357b748930ba Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 08:43:33 -0700 Subject: [PATCH 09/34] add a copy method --- discretize/BaseMesh.py | 6 ++++++ tests/base/test_properties.py | 38 ++++++++++++++++++++++------------- 2 files changed, 30 insertions(+), 14 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index d556ee21a..b5ef59361 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -321,6 +321,12 @@ def save(self, filename='mesh.json', verbose=False): return f + def copy(self): + """ + Make a copy of the current mesh + """ + return properties.copy(self) + class BaseRectangularMesh(BaseMesh): """BaseRectangularMesh""" diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index 8d00027cc..0c2e21c9d 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -3,6 +3,21 @@ import discretize +def compare_meshes(mesh0, mesh1): + assert mesh0.nC == mesh1.nC , ( + 'Number of cells not the same, {} != {}'.format(mesh0.nC, mesh1.nC) + ) + + assert (mesh0.x0 == mesh1.x0).all(), ( + 'x0 different. {} != {}'.format(mesh0.x0, mesh1.x0) + ) + + for i in range(len(mesh0.h)): + assert (mesh0.h[i] == mesh1.h[i]).all(), ( + 'mesh h[{}] different'.format(i) + ) + + class TensorTest(unittest.TestCase): n = [4, 5, 9] @@ -16,21 +31,16 @@ def test_save_load(self): mesh0 = self.mesh f = mesh0.save() mesh1 = discretize.utils.load_mesh(f) - - assert mesh0.nC == mesh1.nC , ( - 'Number of cells not the same, {} != {}'.format(mesh0.nC, mesh1.nC) - ) - - assert (mesh0.x0 == mesh1.x0).all(), ( - 'x0 different. {} != {}'.format(mesh0.x0, mesh1.x0) - ) - - for i in range(len(mesh0.h)): - assert (mesh0.h[i] == mesh1.h[i]).all(), ( - 'mesh h[{}] different'.format(i) - ) - + compare_meshes(mesh0, mesh1) os.remove(f) + print('ok\n') + + def test_copy(self): + print('\nTesting copy of Tensor Mesh ...') + mesh0 = self.mesh + mesh1 = mesh0.copy() + compare_meshes(mesh0, mesh1) + print('ok\n') From 5fc99d5e77fc48fa66b2fbf576460fe4b8e00219 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 09:12:04 -0700 Subject: [PATCH 10/34] fix instantiation of CylMesh so that cartesianOrigin is now a property, add a check to ensure the cartesian origin is being set properly. Have x0 be passed as a kwarg in meshIO --- discretize/BaseMesh.py | 5 +- discretize/CylMesh.py | 39 ++- discretize/MeshIO.py | 2 +- discretize/TreeUtils.c | 613 +++++++++++++++++++++++++++++++++-------- tests/base/test_cyl.py | 6 +- 5 files changed, 537 insertions(+), 128 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index b5ef59361..000e3f782 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -50,7 +50,10 @@ def check_n_shape(self, change): @properties.validator('x0') def check_x0_vs_n(self, change): if len(self.n) != len(change['value']): - raise Exception("Dimension mismatch. x0 != len(n)") + raise Exception( + "Dimension mismatch. x0 has length {} != len(n) which is " + "{}".format(len(x0), len(n)) + ) @property def dim(self): diff --git a/discretize/CylMesh.py b/discretize/CylMesh.py index 05d76c3a8..da647e926 100644 --- a/discretize/CylMesh.py +++ b/discretize/CylMesh.py @@ -1,5 +1,6 @@ from __future__ import print_function import numpy as np +import properties import scipy.sparse as sp from scipy.constants import pi @@ -29,17 +30,39 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): _unitDimensions = [1, 2*np.pi, 1] - def __init__(self, h, x0=None, cartesianOrigin=None): - BaseTensorMesh.__init__(self, h, x0) + cartesianOrigin = properties.Array( + "Cartesian origin of the mesh", + dtype=float, + shape=('*',) + ) + + def __init__(self, h_in=None, **kwargs): + BaseTensorMesh.__init__(self, h_in=h_in, **kwargs) assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" if self.dim == 2: print('Warning, a disk mesh has not been tested thoroughly.') - cartesianOrigin = (np.zeros(self.dim) if cartesianOrigin is None - else cartesianOrigin) - assert len(cartesianOrigin) == self.dim, ("cartesianOrigin must be the " - "same length as the dimension" - " of the mesh.") - self.cartesianOrigin = np.array(cartesianOrigin, dtype=float) + + if 'cartesianOrigin' in kwargs.keys(): + self.cartesianOrigin = kwargs.pop('cartesianOrigin') + else: + self.cartesianOrigin = np.zeros(self.dim) + + # assert len(self.cartesianOrigin) == self.dim, ( + # "cartesianOrigin must be the same length as the dimension" + # " of the mesh." + # ) + # self.cartesianOrigin = np.array(self.cartesianOrigin, dtype=float) + + @properties.validator('cartesianOrigin') + def check_cartesian_origin_shape(self, change): + change['value'] = np.array(change['value'], dtype=float).ravel() + if len(change['value']) != self.dim: + raise Exception( + "Dimension mismatch. The mesh dimension is {}, and the " + "cartesianOrigin provided has length {}".format( + self.dim, len(change['value']) + ) + ) @property def isSymmetric(self): diff --git a/discretize/MeshIO.py b/discretize/MeshIO.py index 42b570609..75a0ef7ea 100644 --- a/discretize/MeshIO.py +++ b/discretize/MeshIO.py @@ -42,7 +42,7 @@ def readCellLine(line): # Adjust the reference point to the bottom south west corner x0[2] = x0[2] - np.sum(h3) # Make the mesh - tensMsh = TensorMesh([h1, h2, h3], x0) + tensMsh = TensorMesh([h1, h2, h3], x0=x0) return tensMsh @classmethod diff --git a/discretize/TreeUtils.c b/discretize/TreeUtils.c index d16088511..e0ae39f81 100644 --- a/discretize/TreeUtils.c +++ b/discretize/TreeUtils.c @@ -1,8 +1,8 @@ -/* Generated by Cython 0.24.1 */ +/* Generated by Cython 0.25.2 */ /* BEGIN: Cython Metadata { - "distutils": {}, + "distutils": {}, "module_name": "discretize.TreeUtils" } END: Cython Metadata */ @@ -14,7 +14,7 @@ END: Cython Metadata */ #elif PY_VERSION_HEX < 0x02060000 || (0x03000000 <= PY_VERSION_HEX && PY_VERSION_HEX < 0x03020000) #error Cython requires Python 2.6+ or Python 3.2+. #else -#define CYTHON_ABI "0_24_1" +#define CYTHON_ABI "0_25_2" #include #ifndef offsetof #define offsetof(type, member) ( (size_t) & ((type*)0) -> member ) @@ -36,6 +36,11 @@ END: Cython Metadata */ #ifndef DL_EXPORT #define DL_EXPORT(t) t #endif +#ifndef HAVE_LONG_LONG + #if PY_VERSION_HEX >= 0x03030000 || (PY_MAJOR_VERSION == 2 && PY_VERSION_HEX >= 0x02070000) + #define HAVE_LONG_LONG + #endif +#endif #ifndef PY_LONG_LONG #define PY_LONG_LONG LONG_LONG #endif @@ -44,13 +49,110 @@ END: Cython Metadata */ #endif #ifdef PYPY_VERSION #define CYTHON_COMPILING_IN_PYPY 1 + #define CYTHON_COMPILING_IN_PYSTON 0 + #define CYTHON_COMPILING_IN_CPYTHON 0 + #undef CYTHON_USE_TYPE_SLOTS + #define CYTHON_USE_TYPE_SLOTS 0 + #undef CYTHON_USE_ASYNC_SLOTS + #define CYTHON_USE_ASYNC_SLOTS 0 + #undef CYTHON_USE_PYLIST_INTERNALS + #define CYTHON_USE_PYLIST_INTERNALS 0 + #undef CYTHON_USE_UNICODE_INTERNALS + #define CYTHON_USE_UNICODE_INTERNALS 0 + #undef CYTHON_USE_UNICODE_WRITER + #define CYTHON_USE_UNICODE_WRITER 0 + #undef CYTHON_USE_PYLONG_INTERNALS + #define CYTHON_USE_PYLONG_INTERNALS 0 + #undef CYTHON_AVOID_BORROWED_REFS + #define CYTHON_AVOID_BORROWED_REFS 1 + #undef CYTHON_ASSUME_SAFE_MACROS + #define CYTHON_ASSUME_SAFE_MACROS 0 + #undef CYTHON_UNPACK_METHODS + #define CYTHON_UNPACK_METHODS 0 + #undef CYTHON_FAST_THREAD_STATE + #define CYTHON_FAST_THREAD_STATE 0 + #undef CYTHON_FAST_PYCALL + #define CYTHON_FAST_PYCALL 0 +#elif defined(PYSTON_VERSION) + #define CYTHON_COMPILING_IN_PYPY 0 + #define CYTHON_COMPILING_IN_PYSTON 1 #define CYTHON_COMPILING_IN_CPYTHON 0 + #ifndef CYTHON_USE_TYPE_SLOTS + #define CYTHON_USE_TYPE_SLOTS 1 + #endif + #undef CYTHON_USE_ASYNC_SLOTS + #define CYTHON_USE_ASYNC_SLOTS 0 + #undef CYTHON_USE_PYLIST_INTERNALS + #define CYTHON_USE_PYLIST_INTERNALS 0 + #ifndef CYTHON_USE_UNICODE_INTERNALS + #define CYTHON_USE_UNICODE_INTERNALS 1 + #endif + #undef CYTHON_USE_UNICODE_WRITER + #define CYTHON_USE_UNICODE_WRITER 0 + #undef CYTHON_USE_PYLONG_INTERNALS + #define CYTHON_USE_PYLONG_INTERNALS 0 + #ifndef CYTHON_AVOID_BORROWED_REFS + #define CYTHON_AVOID_BORROWED_REFS 0 + #endif + #ifndef CYTHON_ASSUME_SAFE_MACROS + #define CYTHON_ASSUME_SAFE_MACROS 1 + #endif + #ifndef CYTHON_UNPACK_METHODS + #define CYTHON_UNPACK_METHODS 1 + #endif + #undef CYTHON_FAST_THREAD_STATE + #define CYTHON_FAST_THREAD_STATE 0 + #undef CYTHON_FAST_PYCALL + #define CYTHON_FAST_PYCALL 0 #else #define CYTHON_COMPILING_IN_PYPY 0 + #define CYTHON_COMPILING_IN_PYSTON 0 #define CYTHON_COMPILING_IN_CPYTHON 1 + #ifndef CYTHON_USE_TYPE_SLOTS + #define CYTHON_USE_TYPE_SLOTS 1 + #endif + #if PY_MAJOR_VERSION < 3 + #undef CYTHON_USE_ASYNC_SLOTS + #define CYTHON_USE_ASYNC_SLOTS 0 + #elif !defined(CYTHON_USE_ASYNC_SLOTS) + #define CYTHON_USE_ASYNC_SLOTS 1 + #endif + #if PY_VERSION_HEX < 0x02070000 + #undef CYTHON_USE_PYLONG_INTERNALS + #define CYTHON_USE_PYLONG_INTERNALS 0 + #elif !defined(CYTHON_USE_PYLONG_INTERNALS) + #define CYTHON_USE_PYLONG_INTERNALS 1 + #endif + #ifndef CYTHON_USE_PYLIST_INTERNALS + #define CYTHON_USE_PYLIST_INTERNALS 1 + #endif + #ifndef CYTHON_USE_UNICODE_INTERNALS + #define CYTHON_USE_UNICODE_INTERNALS 1 + #endif + #if PY_VERSION_HEX < 0x030300F0 + #undef CYTHON_USE_UNICODE_WRITER + #define CYTHON_USE_UNICODE_WRITER 0 + #elif !defined(CYTHON_USE_UNICODE_WRITER) + #define CYTHON_USE_UNICODE_WRITER 1 + #endif + #ifndef CYTHON_AVOID_BORROWED_REFS + #define CYTHON_AVOID_BORROWED_REFS 0 + #endif + #ifndef CYTHON_ASSUME_SAFE_MACROS + #define CYTHON_ASSUME_SAFE_MACROS 1 + #endif + #ifndef CYTHON_UNPACK_METHODS + #define CYTHON_UNPACK_METHODS 1 + #endif + #ifndef CYTHON_FAST_THREAD_STATE + #define CYTHON_FAST_THREAD_STATE 1 + #endif + #ifndef CYTHON_FAST_PYCALL + #define CYTHON_FAST_PYCALL 1 + #endif #endif -#if !defined(CYTHON_USE_PYLONG_INTERNALS) && CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x02070000 - #define CYTHON_USE_PYLONG_INTERNALS 1 +#if !defined(CYTHON_FAST_PYCCALL) +#define CYTHON_FAST_PYCCALL (CYTHON_FAST_PYCALL && PY_VERSION_HEX >= 0x030600B1) #endif #if CYTHON_USE_PYLONG_INTERNALS #include "longintrepr.h" @@ -86,24 +188,44 @@ END: Cython Metadata */ #ifndef Py_TPFLAGS_HAVE_FINALIZE #define Py_TPFLAGS_HAVE_FINALIZE 0 #endif +#ifndef METH_FASTCALL + #define METH_FASTCALL 0x80 + typedef PyObject *(*__Pyx_PyCFunctionFast) (PyObject *self, PyObject **args, + Py_ssize_t nargs, PyObject *kwnames); +#else + #define __Pyx_PyCFunctionFast _PyCFunctionFast +#endif +#if CYTHON_FAST_PYCCALL +#define __Pyx_PyFastCFunction_Check(func)\ + ((PyCFunction_Check(func) && (METH_FASTCALL == (PyCFunction_GET_FLAGS(func) & ~(METH_CLASS | METH_STATIC | METH_COEXIST))))) +#else +#define __Pyx_PyFastCFunction_Check(func) 0 +#endif #if PY_VERSION_HEX > 0x03030000 && defined(PyUnicode_KIND) #define CYTHON_PEP393_ENABLED 1 #define __Pyx_PyUnicode_READY(op) (likely(PyUnicode_IS_READY(op)) ?\ 0 : _PyUnicode_Ready((PyObject *)(op))) #define __Pyx_PyUnicode_GET_LENGTH(u) PyUnicode_GET_LENGTH(u) #define __Pyx_PyUnicode_READ_CHAR(u, i) PyUnicode_READ_CHAR(u, i) + #define __Pyx_PyUnicode_MAX_CHAR_VALUE(u) PyUnicode_MAX_CHAR_VALUE(u) #define __Pyx_PyUnicode_KIND(u) PyUnicode_KIND(u) #define __Pyx_PyUnicode_DATA(u) PyUnicode_DATA(u) #define __Pyx_PyUnicode_READ(k, d, i) PyUnicode_READ(k, d, i) + #define __Pyx_PyUnicode_WRITE(k, d, i, ch) PyUnicode_WRITE(k, d, i, ch) #define __Pyx_PyUnicode_IS_TRUE(u) (0 != (likely(PyUnicode_IS_READY(u)) ? PyUnicode_GET_LENGTH(u) : PyUnicode_GET_SIZE(u))) #else #define CYTHON_PEP393_ENABLED 0 + #define PyUnicode_1BYTE_KIND 1 + #define PyUnicode_2BYTE_KIND 2 + #define PyUnicode_4BYTE_KIND 4 #define __Pyx_PyUnicode_READY(op) (0) #define __Pyx_PyUnicode_GET_LENGTH(u) PyUnicode_GET_SIZE(u) #define __Pyx_PyUnicode_READ_CHAR(u, i) ((Py_UCS4)(PyUnicode_AS_UNICODE(u)[i])) + #define __Pyx_PyUnicode_MAX_CHAR_VALUE(u) ((sizeof(Py_UNICODE) == 2) ? 65535 : 1114111) #define __Pyx_PyUnicode_KIND(u) (sizeof(Py_UNICODE)) #define __Pyx_PyUnicode_DATA(u) ((void*)PyUnicode_AS_UNICODE(u)) #define __Pyx_PyUnicode_READ(k, d, i) ((void)(k), (Py_UCS4)(((Py_UNICODE*)d)[i])) + #define __Pyx_PyUnicode_WRITE(k, d, i, ch) (((void)(k)), ((Py_UNICODE*)d)[i] = ch) #define __Pyx_PyUnicode_IS_TRUE(u) (0 != PyUnicode_GET_SIZE(u)) #endif #if CYTHON_COMPILING_IN_PYPY @@ -128,6 +250,13 @@ END: Cython Metadata */ #define PyObject_Free(p) PyMem_Free(p) #define PyObject_Realloc(p) PyMem_Realloc(p) #endif +#if CYTHON_COMPILING_IN_PYSTON + #define __Pyx_PyCode_HasFreeVars(co) PyCode_HasFreeVars(co) + #define __Pyx_PyFrame_SetLineNumber(frame, lineno) PyFrame_SetLineNumber(frame, lineno) +#else + #define __Pyx_PyCode_HasFreeVars(co) (PyCode_GetNumFree(co) > 0) + #define __Pyx_PyFrame_SetLineNumber(frame, lineno) (frame)->f_lineno = (lineno) +#endif #define __Pyx_PyString_FormatSafe(a, b) ((unlikely((a) == Py_None)) ? PyNumber_Remainder(a, b) : __Pyx_PyString_Format(a, b)) #define __Pyx_PyUnicode_FormatSafe(a, b) ((unlikely((a) == Py_None)) ? PyNumber_Remainder(a, b) : PyUnicode_Format(a, b)) #if PY_MAJOR_VERSION >= 3 @@ -156,6 +285,7 @@ END: Cython Metadata */ #define PySet_CheckExact(obj) (Py_TYPE(obj) == &PySet_Type) #endif #define __Pyx_TypeCheck(obj, type) PyObject_TypeCheck(obj, (PyTypeObject *)type) +#define __Pyx_PyException_Check(obj) __Pyx_TypeCheck(obj, PyExc_Exception) #if PY_MAJOR_VERSION >= 3 #define PyIntObject PyLongObject #define PyInt_Type PyLong_Type @@ -194,18 +324,20 @@ END: Cython Metadata */ #else #define __Pyx_PyMethod_New(func, self, klass) PyMethod_New(func, self, klass) #endif -#if PY_VERSION_HEX >= 0x030500B1 -#define __Pyx_PyAsyncMethodsStruct PyAsyncMethods -#define __Pyx_PyType_AsAsync(obj) (Py_TYPE(obj)->tp_as_async) -#elif CYTHON_COMPILING_IN_CPYTHON && PY_MAJOR_VERSION >= 3 -typedef struct { - unaryfunc am_await; - unaryfunc am_aiter; - unaryfunc am_anext; -} __Pyx_PyAsyncMethodsStruct; -#define __Pyx_PyType_AsAsync(obj) ((__Pyx_PyAsyncMethodsStruct*) (Py_TYPE(obj)->tp_reserved)) +#if CYTHON_USE_ASYNC_SLOTS + #if PY_VERSION_HEX >= 0x030500B1 + #define __Pyx_PyAsyncMethodsStruct PyAsyncMethods + #define __Pyx_PyType_AsAsync(obj) (Py_TYPE(obj)->tp_as_async) + #else + typedef struct { + unaryfunc am_await; + unaryfunc am_aiter; + unaryfunc am_anext; + } __Pyx_PyAsyncMethodsStruct; + #define __Pyx_PyType_AsAsync(obj) ((__Pyx_PyAsyncMethodsStruct*) (Py_TYPE(obj)->tp_reserved)) + #endif #else -#define __Pyx_PyType_AsAsync(obj) NULL + #define __Pyx_PyType_AsAsync(obj) NULL #endif #ifndef CYTHON_RESTRICT #if defined(__GNUC__) @@ -218,10 +350,39 @@ typedef struct { #define CYTHON_RESTRICT #endif #endif +#ifndef CYTHON_UNUSED +# if defined(__GNUC__) +# if !(defined(__cplusplus)) || (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4)) +# define CYTHON_UNUSED __attribute__ ((__unused__)) +# else +# define CYTHON_UNUSED +# endif +# elif defined(__ICC) || (defined(__INTEL_COMPILER) && !defined(_MSC_VER)) +# define CYTHON_UNUSED __attribute__ ((__unused__)) +# else +# define CYTHON_UNUSED +# endif +#endif +#ifndef CYTHON_MAYBE_UNUSED_VAR +# if defined(__cplusplus) + template void CYTHON_MAYBE_UNUSED_VAR( const T& ) { } +# else +# define CYTHON_MAYBE_UNUSED_VAR(x) (void)(x) +# endif +#endif +#ifndef CYTHON_NCP_UNUSED +# if CYTHON_COMPILING_IN_CPYTHON +# define CYTHON_NCP_UNUSED +# else +# define CYTHON_NCP_UNUSED CYTHON_UNUSED +# endif +#endif #define __Pyx_void_to_None(void_result) ((void)(void_result), Py_INCREF(Py_None), Py_None) #ifndef CYTHON_INLINE - #if defined(__GNUC__) + #if defined(__clang__) + #define CYTHON_INLINE __inline__ __attribute__ ((__unused__)) + #elif defined(__GNUC__) #define CYTHON_INLINE __inline__ #elif defined(_MSC_VER) #define CYTHON_INLINE __inline @@ -283,26 +444,6 @@ static CYTHON_INLINE float __PYX_NAN() { #define CYTHON_WITHOUT_ASSERTIONS #endif -#ifndef CYTHON_UNUSED -# if defined(__GNUC__) -# if !(defined(__cplusplus)) || (__GNUC__ > 3 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4)) -# define CYTHON_UNUSED __attribute__ ((__unused__)) -# else -# define CYTHON_UNUSED -# endif -# elif defined(__ICC) || (defined(__INTEL_COMPILER) && !defined(_MSC_VER)) -# define CYTHON_UNUSED __attribute__ ((__unused__)) -# else -# define CYTHON_UNUSED -# endif -#endif -#ifndef CYTHON_NCP_UNUSED -# if CYTHON_COMPILING_IN_CPYTHON -# define CYTHON_NCP_UNUSED -# else -# define CYTHON_NCP_UNUSED CYTHON_UNUSED -# endif -#endif typedef struct {PyObject **p; const char *s; const Py_ssize_t n; const char* encoding; const char is_unicode; const char is_str; const char intern; } __Pyx_StringTabEntry; @@ -380,7 +521,7 @@ static CYTHON_INLINE int __Pyx_PyObject_IsTrue(PyObject*); static CYTHON_INLINE PyObject* __Pyx_PyNumber_IntOrLong(PyObject* x); static CYTHON_INLINE Py_ssize_t __Pyx_PyIndex_AsSsize_t(PyObject*); static CYTHON_INLINE PyObject * __Pyx_PyInt_FromSize_t(size_t); -#if CYTHON_COMPILING_IN_CPYTHON +#if CYTHON_ASSUME_SAFE_MACROS #define __pyx_PyFloat_AsDouble(x) (PyFloat_CheckExact(x) ? PyFloat_AS_DOUBLE(x) : PyFloat_AsDouble(x)) #else #define __pyx_PyFloat_AsDouble(x) PyFloat_AsDouble(x) @@ -560,7 +701,7 @@ static const char *__pyx_f[] = { #define __Pyx_XCLEAR(r) do { if((r) != NULL) {PyObject* tmp = ((PyObject*)(r)); r = NULL; __Pyx_DECREF(tmp);}} while(0) /* PyObjectGetAttrStr.proto */ -#if CYTHON_COMPILING_IN_CPYTHON +#if CYTHON_USE_TYPE_SLOTS static CYTHON_INLINE PyObject* __Pyx_PyObject_GetAttrStr(PyObject* obj, PyObject* attr_name) { PyTypeObject* tp = Py_TYPE(obj); if (likely(tp->tp_getattro)) @@ -595,7 +736,7 @@ static CYTHON_INLINE int __Pyx_ArgTypeTest(PyObject *obj, PyTypeObject *type, in const char *name, int exact); /* ListCompAppend.proto */ -#if CYTHON_COMPILING_IN_CPYTHON +#if CYTHON_USE_PYLIST_INTERNALS && CYTHON_ASSUME_SAFE_MACROS static CYTHON_INLINE int __Pyx_ListComp_Append(PyObject* list, PyObject* x) { PyListObject* L = (PyListObject*) list; Py_ssize_t len = Py_SIZE(list); @@ -646,6 +787,24 @@ static CYTHON_INLINE PyObject *__Pyx_GetItemInt_Generic(PyObject *o, PyObject* j static CYTHON_INLINE PyObject *__Pyx_GetItemInt_Fast(PyObject *o, Py_ssize_t i, int is_list, int wraparound, int boundscheck); +/* PyFunctionFastCall.proto */ +#if CYTHON_FAST_PYCALL +#define __Pyx_PyFunction_FastCall(func, args, nargs)\ + __Pyx_PyFunction_FastCallDict((func), (args), (nargs), NULL) +#if 1 || PY_VERSION_HEX < 0x030600B1 +static PyObject *__Pyx_PyFunction_FastCallDict(PyObject *func, PyObject **args, int nargs, PyObject *kwargs); +#else +#define __Pyx_PyFunction_FastCallDict(func, args, nargs, kwargs) _PyFunction_FastCallDict(func, args, nargs, kwargs) +#endif +#endif + +/* PyCFunctionFastCall.proto */ +#if CYTHON_FAST_PYCCALL +static CYTHON_INLINE PyObject *__Pyx_PyCFunction_FastCall(PyObject *func, PyObject **args, Py_ssize_t nargs); +#else +#define __Pyx_PyCFunction_FastCall(func, args, nargs) (assert(0), NULL) +#endif + /* PyObjectCall.proto */ #if CYTHON_COMPILING_IN_CPYTHON static CYTHON_INLINE PyObject* __Pyx_PyObject_Call(PyObject *func, PyObject *arg, PyObject *kw); @@ -744,8 +903,8 @@ static const char __pyx_k_dimension[] = "dimension"; static const char __pyx_k_levelBits[] = "levelBits"; static const char __pyx_k_discretize_TreeUtils[] = "discretize.TreeUtils"; static const char __pyx_k_The_Z_order_curve_is_generated[] = "\n The Z-order curve is generated by interleaving the bits of an offset.\n\n See also:\n\n https://github.com/cortesi/scurve\n Aldo Cortesi \n\n"; -static const char __pyx_k_Users_josephcapriotti_python_pa[] = "/Users/josephcapriotti/python_packages/discretize/discretize/TreeUtils.pyx"; -static PyObject *__pyx_kp_s_Users_josephcapriotti_python_pa; +static const char __pyx_k_Users_lindseyjh_git_simpeg_disc[] = "/Users/lindseyjh/git/simpeg/discretize/discretize/TreeUtils.pyx"; +static PyObject *__pyx_kp_s_Users_lindseyjh_git_simpeg_disc; static PyObject *__pyx_n_s__3; static PyObject *__pyx_n_s_b; static PyObject *__pyx_n_s_bitoff; @@ -1026,8 +1185,9 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_2index(CYTHON_UNUSED PyObject PyObject *__pyx_t_9 = NULL; PyObject *__pyx_t_10 = NULL; PyObject *__pyx_t_11 = NULL; - PyObject *__pyx_t_12 = NULL; - PY_LONG_LONG __pyx_t_13; + int __pyx_t_12; + PyObject *__pyx_t_13 = NULL; + PY_LONG_LONG __pyx_t_14; __Pyx_RefNannySetupContext("index", 0); __Pyx_INCREF(__pyx_v_p); @@ -1056,7 +1216,7 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_2index(CYTHON_UNUSED PyObject __pyx_t_2 = __pyx_v_p; __Pyx_INCREF(__pyx_t_2); __pyx_t_3 = 0; for (;;) { if (__pyx_t_3 >= PyList_GET_SIZE(__pyx_t_2)) break; - #if CYTHON_COMPILING_IN_CPYTHON + #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS __pyx_t_4 = PyList_GET_ITEM(__pyx_t_2, __pyx_t_3); __Pyx_INCREF(__pyx_t_4); __pyx_t_3++; if (unlikely(0 < 0)) __PYX_ERR(0, 33, __pyx_L1_error) #else __pyx_t_4 = PySequence_ITEM(__pyx_t_2, __pyx_t_3); __pyx_t_3++; if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 33, __pyx_L1_error) @@ -1147,47 +1307,73 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_2index(CYTHON_UNUSED PyObject __pyx_t_10 = __Pyx_PyInt_From_long((__pyx_v_bitoff + 1)); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 40, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_t_11 = NULL; - __pyx_t_3 = 0; - if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_2))) { + __pyx_t_12 = 0; + if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_2))) { __pyx_t_11 = PyMethod_GET_SELF(__pyx_t_2); if (likely(__pyx_t_11)) { PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2); __Pyx_INCREF(__pyx_t_11); __Pyx_INCREF(function); __Pyx_DECREF_SET(__pyx_t_2, function); - __pyx_t_3 = 1; + __pyx_t_12 = 1; } } - __pyx_t_12 = PyTuple_New(4+__pyx_t_3); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 40, __pyx_L1_error) - __Pyx_GOTREF(__pyx_t_12); - if (__pyx_t_11) { - __Pyx_GIVEREF(__pyx_t_11); PyTuple_SET_ITEM(__pyx_t_12, 0, __pyx_t_11); __pyx_t_11 = NULL; - } - __Pyx_GIVEREF(__pyx_t_4); - PyTuple_SET_ITEM(__pyx_t_12, 0+__pyx_t_3, __pyx_t_4); - __Pyx_GIVEREF(__pyx_t_8); - PyTuple_SET_ITEM(__pyx_t_12, 1+__pyx_t_3, __pyx_t_8); - __Pyx_GIVEREF(__pyx_t_9); - PyTuple_SET_ITEM(__pyx_t_12, 2+__pyx_t_3, __pyx_t_9); - __Pyx_GIVEREF(__pyx_t_10); - PyTuple_SET_ITEM(__pyx_t_12, 3+__pyx_t_3, __pyx_t_10); - __pyx_t_4 = 0; - __pyx_t_8 = 0; - __pyx_t_9 = 0; - __pyx_t_10 = 0; - __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_12, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 40, __pyx_L1_error) - __Pyx_GOTREF(__pyx_t_1); - __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0; + #if CYTHON_FAST_PYCALL + if (PyFunction_Check(__pyx_t_2)) { + PyObject *__pyx_temp[5] = {__pyx_t_11, __pyx_t_4, __pyx_t_8, __pyx_t_9, __pyx_t_10}; + __pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_2, __pyx_temp+1-__pyx_t_12, 4+__pyx_t_12); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 40, __pyx_L1_error) + __Pyx_XDECREF(__pyx_t_11); __pyx_t_11 = 0; + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; + __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0; + __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; + __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0; + } else + #endif + #if CYTHON_FAST_PYCCALL + if (__Pyx_PyFastCFunction_Check(__pyx_t_2)) { + PyObject *__pyx_temp[5] = {__pyx_t_11, __pyx_t_4, __pyx_t_8, __pyx_t_9, __pyx_t_10}; + __pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_2, __pyx_temp+1-__pyx_t_12, 4+__pyx_t_12); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 40, __pyx_L1_error) + __Pyx_XDECREF(__pyx_t_11); __pyx_t_11 = 0; + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; + __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0; + __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; + __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0; + } else + #endif + { + __pyx_t_13 = PyTuple_New(4+__pyx_t_12); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 40, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_13); + if (__pyx_t_11) { + __Pyx_GIVEREF(__pyx_t_11); PyTuple_SET_ITEM(__pyx_t_13, 0, __pyx_t_11); __pyx_t_11 = NULL; + } + __Pyx_GIVEREF(__pyx_t_4); + PyTuple_SET_ITEM(__pyx_t_13, 0+__pyx_t_12, __pyx_t_4); + __Pyx_GIVEREF(__pyx_t_8); + PyTuple_SET_ITEM(__pyx_t_13, 1+__pyx_t_12, __pyx_t_8); + __Pyx_GIVEREF(__pyx_t_9); + PyTuple_SET_ITEM(__pyx_t_13, 2+__pyx_t_12, __pyx_t_9); + __Pyx_GIVEREF(__pyx_t_10); + PyTuple_SET_ITEM(__pyx_t_13, 3+__pyx_t_12, __pyx_t_10); + __pyx_t_4 = 0; + __pyx_t_8 = 0; + __pyx_t_9 = 0; + __pyx_t_10 = 0; + __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_13, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 40, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_13); __pyx_t_13 = 0; + } __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; __pyx_t_2 = __Pyx_PyInt_From_int(__pyx_v_i); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 40, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_2); - __pyx_t_12 = PyNumber_Lshift(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 40, __pyx_L1_error) - __Pyx_GOTREF(__pyx_t_12); + __pyx_t_13 = PyNumber_Lshift(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 40, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_13); __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; - __pyx_t_13 = __Pyx_PyInt_As_PY_LONG_LONG(__pyx_t_12); if (unlikely((__pyx_t_13 == (PY_LONG_LONG)-1) && PyErr_Occurred())) __PYX_ERR(0, 40, __pyx_L1_error) - __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0; - __pyx_v_b = __pyx_t_13; + __pyx_t_14 = __Pyx_PyInt_As_PY_LONG_LONG(__pyx_t_13); if (unlikely((__pyx_t_14 == (PY_LONG_LONG)-1) && PyErr_Occurred())) __PYX_ERR(0, 40, __pyx_L1_error) + __Pyx_DECREF(__pyx_t_13); __pyx_t_13 = 0; + __pyx_v_b = __pyx_t_14; /* "discretize/TreeUtils.pyx":41 * poff = dimension-(i % dimension)-1 @@ -1207,10 +1393,10 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_2index(CYTHON_UNUSED PyObject * */ __Pyx_XDECREF(__pyx_r); - __pyx_t_12 = __Pyx_PyInt_From_PY_LONG_LONG(((__pyx_v_idx << __pyx_v_levelBits) + __pyx_v_level)); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 43, __pyx_L1_error) - __Pyx_GOTREF(__pyx_t_12); - __pyx_r = __pyx_t_12; - __pyx_t_12 = 0; + __pyx_t_13 = __Pyx_PyInt_From_PY_LONG_LONG(((__pyx_v_idx << __pyx_v_levelBits) + __pyx_v_level)); if (unlikely(!__pyx_t_13)) __PYX_ERR(0, 43, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_13); + __pyx_r = __pyx_t_13; + __pyx_t_13 = 0; goto __pyx_L0; /* "discretize/TreeUtils.pyx":26 @@ -1230,7 +1416,7 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_2index(CYTHON_UNUSED PyObject __Pyx_XDECREF(__pyx_t_9); __Pyx_XDECREF(__pyx_t_10); __Pyx_XDECREF(__pyx_t_11); - __Pyx_XDECREF(__pyx_t_12); + __Pyx_XDECREF(__pyx_t_13); __Pyx_AddTraceback("discretize.TreeUtils.index", __pyx_clineno, __pyx_lineno, __pyx_filename); __pyx_r = NULL; __pyx_L0:; @@ -1343,12 +1529,11 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_4point(CYTHON_UNUSED PyObject PyObject *__pyx_t_7 = NULL; PyObject *__pyx_t_8 = NULL; PyObject *__pyx_t_9 = NULL; - Py_ssize_t __pyx_t_10; + int __pyx_t_10; PyObject *__pyx_t_11 = NULL; long __pyx_t_12; PY_LONG_LONG __pyx_t_13; int __pyx_t_14; - int __pyx_t_15; __Pyx_RefNannySetupContext("point", 0); /* "discretize/TreeUtils.pyx":52 @@ -1427,7 +1612,7 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_4point(CYTHON_UNUSED PyObject __Pyx_GOTREF(__pyx_t_8); __pyx_t_9 = NULL; __pyx_t_10 = 0; - if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_4))) { + if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) { __pyx_t_9 = PyMethod_GET_SELF(__pyx_t_4); if (likely(__pyx_t_9)) { PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4); @@ -1437,26 +1622,52 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_4point(CYTHON_UNUSED PyObject __pyx_t_10 = 1; } } - __pyx_t_11 = PyTuple_New(4+__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 58, __pyx_L1_error) - __Pyx_GOTREF(__pyx_t_11); - if (__pyx_t_9) { - __Pyx_GIVEREF(__pyx_t_9); PyTuple_SET_ITEM(__pyx_t_11, 0, __pyx_t_9); __pyx_t_9 = NULL; - } - __Pyx_GIVEREF(__pyx_t_5); - PyTuple_SET_ITEM(__pyx_t_11, 0+__pyx_t_10, __pyx_t_5); - __Pyx_GIVEREF(__pyx_t_6); - PyTuple_SET_ITEM(__pyx_t_11, 1+__pyx_t_10, __pyx_t_6); - __Pyx_GIVEREF(__pyx_t_7); - PyTuple_SET_ITEM(__pyx_t_11, 2+__pyx_t_10, __pyx_t_7); - __Pyx_GIVEREF(__pyx_t_8); - PyTuple_SET_ITEM(__pyx_t_11, 3+__pyx_t_10, __pyx_t_8); - __pyx_t_5 = 0; - __pyx_t_6 = 0; - __pyx_t_7 = 0; - __pyx_t_8 = 0; - __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_11, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 58, __pyx_L1_error) - __Pyx_GOTREF(__pyx_t_1); - __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0; + #if CYTHON_FAST_PYCALL + if (PyFunction_Check(__pyx_t_4)) { + PyObject *__pyx_temp[5] = {__pyx_t_9, __pyx_t_5, __pyx_t_6, __pyx_t_7, __pyx_t_8}; + __pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-__pyx_t_10, 4+__pyx_t_10); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 58, __pyx_L1_error) + __Pyx_XDECREF(__pyx_t_9); __pyx_t_9 = 0; + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; + __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; + __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0; + __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0; + } else + #endif + #if CYTHON_FAST_PYCCALL + if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) { + PyObject *__pyx_temp[5] = {__pyx_t_9, __pyx_t_5, __pyx_t_6, __pyx_t_7, __pyx_t_8}; + __pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-__pyx_t_10, 4+__pyx_t_10); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 58, __pyx_L1_error) + __Pyx_XDECREF(__pyx_t_9); __pyx_t_9 = 0; + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; + __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; + __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0; + __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0; + } else + #endif + { + __pyx_t_11 = PyTuple_New(4+__pyx_t_10); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 58, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_11); + if (__pyx_t_9) { + __Pyx_GIVEREF(__pyx_t_9); PyTuple_SET_ITEM(__pyx_t_11, 0, __pyx_t_9); __pyx_t_9 = NULL; + } + __Pyx_GIVEREF(__pyx_t_5); + PyTuple_SET_ITEM(__pyx_t_11, 0+__pyx_t_10, __pyx_t_5); + __Pyx_GIVEREF(__pyx_t_6); + PyTuple_SET_ITEM(__pyx_t_11, 1+__pyx_t_10, __pyx_t_6); + __Pyx_GIVEREF(__pyx_t_7); + PyTuple_SET_ITEM(__pyx_t_11, 2+__pyx_t_10, __pyx_t_7); + __Pyx_GIVEREF(__pyx_t_8); + PyTuple_SET_ITEM(__pyx_t_11, 3+__pyx_t_10, __pyx_t_8); + __pyx_t_5 = 0; + __pyx_t_6 = 0; + __pyx_t_7 = 0; + __pyx_t_8 = 0; + __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_11, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 58, __pyx_L1_error) + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0; + } __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __pyx_t_12 = ((__pyx_v_iwidth - __pyx_v_i) - 1); if (unlikely(__pyx_v_dimension == 0)) { @@ -1488,8 +1699,8 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_4point(CYTHON_UNUSED PyObject PyErr_SetString(PyExc_ZeroDivisionError, "integer division or modulo by zero"); __PYX_ERR(0, 59, __pyx_L1_error) } - __pyx_t_14 = __Pyx_mod_int(__pyx_v_i, __pyx_v_dimension); - __pyx_t_11 = __Pyx_GetItemInt_List(__pyx_v_p, __pyx_t_14, int, 1, __Pyx_PyInt_From_int, 1, 1, 1); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 59, __pyx_L1_error) + __pyx_t_10 = __Pyx_mod_int(__pyx_v_i, __pyx_v_dimension); + __pyx_t_11 = __Pyx_GetItemInt_List(__pyx_v_p, __pyx_t_10, int, 1, __Pyx_PyInt_From_int, 1, 1, 1); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 59, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __pyx_t_4 = __Pyx_PyInt_From_PY_LONG_LONG(__pyx_v_b); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 59, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); @@ -1497,7 +1708,7 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_4point(CYTHON_UNUSED PyObject __Pyx_GOTREF(__pyx_t_1); __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; - if (unlikely(__Pyx_SetItemInt(__pyx_v_p, __pyx_t_14, __pyx_t_1, int, 1, __Pyx_PyInt_From_int, 1, 1, 1) < 0)) __PYX_ERR(0, 59, __pyx_L1_error) + if (unlikely(__Pyx_SetItemInt(__pyx_v_p, __pyx_t_10, __pyx_t_1, int, 1, __Pyx_PyInt_From_int, 1, 1, 1) < 0)) __PYX_ERR(0, 59, __pyx_L1_error) __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; } @@ -1508,7 +1719,7 @@ static PyObject *__pyx_pf_10discretize_9TreeUtils_4point(CYTHON_UNUSED PyObject * return p + [n] * */ - __pyx_t_15 = PyList_Reverse(__pyx_v_p); if (unlikely(__pyx_t_15 == -1)) __PYX_ERR(0, 60, __pyx_L1_error) + __pyx_t_14 = PyList_Reverse(__pyx_v_p); if (unlikely(__pyx_t_14 == -1)) __PYX_ERR(0, 60, __pyx_L1_error) /* "discretize/TreeUtils.pyx":61 * p[i % dimension] |= b @@ -1582,7 +1793,7 @@ static struct PyModuleDef __pyx_moduledef = { #endif static __Pyx_StringTabEntry __pyx_string_tab[] = { - {&__pyx_kp_s_Users_josephcapriotti_python_pa, __pyx_k_Users_josephcapriotti_python_pa, sizeof(__pyx_k_Users_josephcapriotti_python_pa), 0, 0, 1, 0}, + {&__pyx_kp_s_Users_lindseyjh_git_simpeg_disc, __pyx_k_Users_lindseyjh_git_simpeg_disc, sizeof(__pyx_k_Users_lindseyjh_git_simpeg_disc), 0, 0, 1, 0}, {&__pyx_n_s__3, __pyx_k__3, sizeof(__pyx_k__3), 0, 0, 1, 1}, {&__pyx_n_s_b, __pyx_k_b, sizeof(__pyx_k_b), 0, 0, 1, 1}, {&__pyx_n_s_bitoff, __pyx_k_bitoff, sizeof(__pyx_k_bitoff), 0, 0, 1, 1}, @@ -1630,7 +1841,7 @@ static int __Pyx_InitCachedConstants(void) { __pyx_tuple_ = PyTuple_Pack(4, __pyx_n_s_x, __pyx_n_s_width, __pyx_n_s_start, __pyx_n_s_end); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 18, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple_); __Pyx_GIVEREF(__pyx_tuple_); - __pyx_codeobj__2 = (PyObject*)__Pyx_PyCode_New(4, 0, 4, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple_, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_josephcapriotti_python_pa, __pyx_n_s_bitrange, 18, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__2)) __PYX_ERR(0, 18, __pyx_L1_error) + __pyx_codeobj__2 = (PyObject*)__Pyx_PyCode_New(4, 0, 4, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple_, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_lindseyjh_git_simpeg_disc, __pyx_n_s_bitrange, 18, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__2)) __PYX_ERR(0, 18, __pyx_L1_error) /* "discretize/TreeUtils.pyx":26 * @@ -1642,7 +1853,7 @@ static int __Pyx_InitCachedConstants(void) { __pyx_tuple__4 = PyTuple_Pack(12, __pyx_n_s_dimension, __pyx_n_s_bits, __pyx_n_s_levelBits, __pyx_n_s_p, __pyx_n_s_level, __pyx_n_s_idx, __pyx_n_s_iwidth, __pyx_n_s_i, __pyx_n_s_b, __pyx_n_s_bitoff, __pyx_n_s_poff, __pyx_n_s__3); if (unlikely(!__pyx_tuple__4)) __PYX_ERR(0, 26, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__4); __Pyx_GIVEREF(__pyx_tuple__4); - __pyx_codeobj__5 = (PyObject*)__Pyx_PyCode_New(5, 0, 12, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__4, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_josephcapriotti_python_pa, __pyx_n_s_index, 26, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__5)) __PYX_ERR(0, 26, __pyx_L1_error) + __pyx_codeobj__5 = (PyObject*)__Pyx_PyCode_New(5, 0, 12, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__4, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_lindseyjh_git_simpeg_disc, __pyx_n_s_index, 26, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__5)) __PYX_ERR(0, 26, __pyx_L1_error) /* "discretize/TreeUtils.pyx":46 * @@ -1654,7 +1865,7 @@ static int __Pyx_InitCachedConstants(void) { __pyx_tuple__6 = PyTuple_Pack(9, __pyx_n_s_dimension, __pyx_n_s_bits, __pyx_n_s_levelBits, __pyx_n_s_idx, __pyx_n_s_p, __pyx_n_s_iwidth, __pyx_n_s_i, __pyx_n_s_n, __pyx_n_s_b); if (unlikely(!__pyx_tuple__6)) __PYX_ERR(0, 46, __pyx_L1_error) __Pyx_GOTREF(__pyx_tuple__6); __Pyx_GIVEREF(__pyx_tuple__6); - __pyx_codeobj__7 = (PyObject*)__Pyx_PyCode_New(4, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__6, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_josephcapriotti_python_pa, __pyx_n_s_point, 46, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__7)) __PYX_ERR(0, 46, __pyx_L1_error) + __pyx_codeobj__7 = (PyObject*)__Pyx_PyCode_New(4, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__6, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_lindseyjh_git_simpeg_disc, __pyx_n_s_point, 46, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__7)) __PYX_ERR(0, 46, __pyx_L1_error) __Pyx_RefNannyFinishContext(); return 0; __pyx_L1_error:; @@ -2049,7 +2260,7 @@ static CYTHON_INLINE int __Pyx_mod_int(int a, int b) { /* GetModuleGlobalName */ static CYTHON_INLINE PyObject *__Pyx_GetModuleGlobalName(PyObject *name) { PyObject *result; -#if CYTHON_COMPILING_IN_CPYTHON +#if !CYTHON_AVOID_BORROWED_REFS result = PyDict_GetItem(__pyx_d, name); if (likely(result)) { Py_INCREF(result); @@ -2075,7 +2286,7 @@ static CYTHON_INLINE PyObject *__Pyx_GetModuleGlobalName(PyObject *name) { static CYTHON_INLINE PyObject *__Pyx_GetItemInt_List_Fast(PyObject *o, Py_ssize_t i, CYTHON_NCP_UNUSED int wraparound, CYTHON_NCP_UNUSED int boundscheck) { -#if CYTHON_COMPILING_IN_CPYTHON +#if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS if (wraparound & unlikely(i < 0)) i += PyList_GET_SIZE(o); if ((!boundscheck) || likely((0 <= i) & (i < PyList_GET_SIZE(o)))) { PyObject *r = PyList_GET_ITEM(o, i); @@ -2090,7 +2301,7 @@ static CYTHON_INLINE PyObject *__Pyx_GetItemInt_List_Fast(PyObject *o, Py_ssize_ static CYTHON_INLINE PyObject *__Pyx_GetItemInt_Tuple_Fast(PyObject *o, Py_ssize_t i, CYTHON_NCP_UNUSED int wraparound, CYTHON_NCP_UNUSED int boundscheck) { -#if CYTHON_COMPILING_IN_CPYTHON +#if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS if (wraparound & unlikely(i < 0)) i += PyTuple_GET_SIZE(o); if ((!boundscheck) || likely((0 <= i) & (i < PyTuple_GET_SIZE(o)))) { PyObject *r = PyTuple_GET_ITEM(o, i); @@ -2105,7 +2316,7 @@ static CYTHON_INLINE PyObject *__Pyx_GetItemInt_Tuple_Fast(PyObject *o, Py_ssize static CYTHON_INLINE PyObject *__Pyx_GetItemInt_Fast(PyObject *o, Py_ssize_t i, int is_list, CYTHON_NCP_UNUSED int wraparound, CYTHON_NCP_UNUSED int boundscheck) { -#if CYTHON_COMPILING_IN_CPYTHON +#if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS && CYTHON_USE_TYPE_SLOTS if (is_list || PyList_CheckExact(o)) { Py_ssize_t n = ((!wraparound) | likely(i >= 0)) ? i : i + PyList_GET_SIZE(o); if ((!boundscheck) || (likely((n >= 0) & (n < PyList_GET_SIZE(o))))) { @@ -2145,6 +2356,144 @@ static CYTHON_INLINE PyObject *__Pyx_GetItemInt_Fast(PyObject *o, Py_ssize_t i, return __Pyx_GetItemInt_Generic(o, PyInt_FromSsize_t(i)); } +/* PyFunctionFastCall */ + #if CYTHON_FAST_PYCALL +#include "frameobject.h" +static PyObject* __Pyx_PyFunction_FastCallNoKw(PyCodeObject *co, PyObject **args, Py_ssize_t na, + PyObject *globals) { + PyFrameObject *f; + PyThreadState *tstate = PyThreadState_GET(); + PyObject **fastlocals; + Py_ssize_t i; + PyObject *result; + assert(globals != NULL); + /* XXX Perhaps we should create a specialized + PyFrame_New() that doesn't take locals, but does + take builtins without sanity checking them. + */ + assert(tstate != NULL); + f = PyFrame_New(tstate, co, globals, NULL); + if (f == NULL) { + return NULL; + } + fastlocals = f->f_localsplus; + for (i = 0; i < na; i++) { + Py_INCREF(*args); + fastlocals[i] = *args++; + } + result = PyEval_EvalFrameEx(f,0); + ++tstate->recursion_depth; + Py_DECREF(f); + --tstate->recursion_depth; + return result; +} +#if 1 || PY_VERSION_HEX < 0x030600B1 +static PyObject *__Pyx_PyFunction_FastCallDict(PyObject *func, PyObject **args, int nargs, PyObject *kwargs) { + PyCodeObject *co = (PyCodeObject *)PyFunction_GET_CODE(func); + PyObject *globals = PyFunction_GET_GLOBALS(func); + PyObject *argdefs = PyFunction_GET_DEFAULTS(func); + PyObject *closure; +#if PY_MAJOR_VERSION >= 3 + PyObject *kwdefs; +#endif + PyObject *kwtuple, **k; + PyObject **d; + Py_ssize_t nd; + Py_ssize_t nk; + PyObject *result; + assert(kwargs == NULL || PyDict_Check(kwargs)); + nk = kwargs ? PyDict_Size(kwargs) : 0; + if (Py_EnterRecursiveCall((char*)" while calling a Python object")) { + return NULL; + } + if ( +#if PY_MAJOR_VERSION >= 3 + co->co_kwonlyargcount == 0 && +#endif + likely(kwargs == NULL || nk == 0) && + co->co_flags == (CO_OPTIMIZED | CO_NEWLOCALS | CO_NOFREE)) { + if (argdefs == NULL && co->co_argcount == nargs) { + result = __Pyx_PyFunction_FastCallNoKw(co, args, nargs, globals); + goto done; + } + else if (nargs == 0 && argdefs != NULL + && co->co_argcount == Py_SIZE(argdefs)) { + /* function called with no arguments, but all parameters have + a default value: use default values as arguments .*/ + args = &PyTuple_GET_ITEM(argdefs, 0); + result =__Pyx_PyFunction_FastCallNoKw(co, args, Py_SIZE(argdefs), globals); + goto done; + } + } + if (kwargs != NULL) { + Py_ssize_t pos, i; + kwtuple = PyTuple_New(2 * nk); + if (kwtuple == NULL) { + result = NULL; + goto done; + } + k = &PyTuple_GET_ITEM(kwtuple, 0); + pos = i = 0; + while (PyDict_Next(kwargs, &pos, &k[i], &k[i+1])) { + Py_INCREF(k[i]); + Py_INCREF(k[i+1]); + i += 2; + } + nk = i / 2; + } + else { + kwtuple = NULL; + k = NULL; + } + closure = PyFunction_GET_CLOSURE(func); +#if PY_MAJOR_VERSION >= 3 + kwdefs = PyFunction_GET_KW_DEFAULTS(func); +#endif + if (argdefs != NULL) { + d = &PyTuple_GET_ITEM(argdefs, 0); + nd = Py_SIZE(argdefs); + } + else { + d = NULL; + nd = 0; + } +#if PY_MAJOR_VERSION >= 3 + result = PyEval_EvalCodeEx((PyObject*)co, globals, (PyObject *)NULL, + args, nargs, + k, (int)nk, + d, (int)nd, kwdefs, closure); +#else + result = PyEval_EvalCodeEx(co, globals, (PyObject *)NULL, + args, nargs, + k, (int)nk, + d, (int)nd, closure); +#endif + Py_XDECREF(kwtuple); +done: + Py_LeaveRecursiveCall(); + return result; +} +#endif // CPython < 3.6 +#endif // CYTHON_FAST_PYCALL + +/* PyCFunctionFastCall */ + #if CYTHON_FAST_PYCCALL +static CYTHON_INLINE PyObject * __Pyx_PyCFunction_FastCall(PyObject *func_obj, PyObject **args, Py_ssize_t nargs) { + PyCFunctionObject *func = (PyCFunctionObject*)func_obj; + PyCFunction meth = PyCFunction_GET_FUNCTION(func); + PyObject *self = PyCFunction_GET_SELF(func); + assert(PyCFunction_Check(func)); + assert(METH_FASTCALL == (PyCFunction_GET_FLAGS(func) & ~(METH_CLASS | METH_STATIC | METH_COEXIST))); + assert(nargs >= 0); + assert(nargs == 0 || args != NULL); + /* _PyCFunction_FastCallDict() must not be called with an exception set, + because it may clear it (directly or indirectly) and so the + caller loses its exception */ + assert(!PyErr_Occurred()); + return (*((__Pyx_PyCFunctionFast)meth)) (self, args, nargs, NULL); +} +#endif // CYTHON_FAST_PYCCALL + /* PyObjectCall */ #if CYTHON_COMPILING_IN_CPYTHON static CYTHON_INLINE PyObject* __Pyx_PyObject_Call(PyObject *func, PyObject *arg, PyObject *kw) { @@ -2183,7 +2532,7 @@ static CYTHON_INLINE PyObject* __Pyx_PyObject_Call(PyObject *func, PyObject *arg } static CYTHON_INLINE int __Pyx_SetItemInt_Fast(PyObject *o, Py_ssize_t i, PyObject *v, int is_list, CYTHON_NCP_UNUSED int wraparound, CYTHON_NCP_UNUSED int boundscheck) { -#if CYTHON_COMPILING_IN_CPYTHON +#if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS && CYTHON_USE_TYPE_SLOTS if (is_list || PyList_CheckExact(o)) { Py_ssize_t n = (!wraparound) ? i : ((likely(i >= 0)) ? i : i + PyList_GET_SIZE(o)); if ((!boundscheck) || likely((n >= 0) & (n < PyList_GET_SIZE(o)))) { @@ -2375,7 +2724,7 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, 0 /*PyObject *locals*/ ); if (!py_frame) goto bad; - py_frame->f_lineno = py_line; + __Pyx_PyFrame_SetLineNumber(py_frame, py_line); PyTraceBack_Here(py_frame); bad: Py_XDECREF(py_code); @@ -2438,14 +2787,18 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, return PyInt_FromLong((long) value); } else if (sizeof(PY_LONG_LONG) <= sizeof(unsigned long)) { return PyLong_FromUnsignedLong((unsigned long) value); +#ifdef HAVE_LONG_LONG } else if (sizeof(PY_LONG_LONG) <= sizeof(unsigned PY_LONG_LONG)) { return PyLong_FromUnsignedLongLong((unsigned PY_LONG_LONG) value); +#endif } } else { if (sizeof(PY_LONG_LONG) <= sizeof(long)) { return PyInt_FromLong((long) value); +#ifdef HAVE_LONG_LONG } else if (sizeof(PY_LONG_LONG) <= sizeof(PY_LONG_LONG)) { return PyLong_FromLongLong((PY_LONG_LONG) value); +#endif } } { @@ -2465,14 +2818,18 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, return PyInt_FromLong((long) value); } else if (sizeof(int) <= sizeof(unsigned long)) { return PyLong_FromUnsignedLong((unsigned long) value); +#ifdef HAVE_LONG_LONG } else if (sizeof(int) <= sizeof(unsigned PY_LONG_LONG)) { return PyLong_FromUnsignedLongLong((unsigned PY_LONG_LONG) value); +#endif } } else { if (sizeof(int) <= sizeof(long)) { return PyInt_FromLong((long) value); +#ifdef HAVE_LONG_LONG } else if (sizeof(int) <= sizeof(PY_LONG_LONG)) { return PyLong_FromLongLong((PY_LONG_LONG) value); +#endif } } { @@ -2492,14 +2849,18 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, return PyInt_FromLong((long) value); } else if (sizeof(long) <= sizeof(unsigned long)) { return PyLong_FromUnsignedLong((unsigned long) value); +#ifdef HAVE_LONG_LONG } else if (sizeof(long) <= sizeof(unsigned PY_LONG_LONG)) { return PyLong_FromUnsignedLongLong((unsigned PY_LONG_LONG) value); +#endif } } else { if (sizeof(long) <= sizeof(long)) { return PyInt_FromLong((long) value); +#ifdef HAVE_LONG_LONG } else if (sizeof(long) <= sizeof(PY_LONG_LONG)) { return PyLong_FromLongLong((PY_LONG_LONG) value); +#endif } } { @@ -2578,8 +2939,10 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, #endif if (sizeof(PY_LONG_LONG) <= sizeof(unsigned long)) { __PYX_VERIFY_RETURN_INT_EXC(PY_LONG_LONG, unsigned long, PyLong_AsUnsignedLong(x)) +#ifdef HAVE_LONG_LONG } else if (sizeof(PY_LONG_LONG) <= sizeof(unsigned PY_LONG_LONG)) { __PYX_VERIFY_RETURN_INT_EXC(PY_LONG_LONG, unsigned PY_LONG_LONG, PyLong_AsUnsignedLongLong(x)) +#endif } } else { #if CYTHON_USE_PYLONG_INTERNALS @@ -2646,8 +3009,10 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, #endif if (sizeof(PY_LONG_LONG) <= sizeof(long)) { __PYX_VERIFY_RETURN_INT_EXC(PY_LONG_LONG, long, PyLong_AsLong(x)) +#ifdef HAVE_LONG_LONG } else if (sizeof(PY_LONG_LONG) <= sizeof(PY_LONG_LONG)) { __PYX_VERIFY_RETURN_INT_EXC(PY_LONG_LONG, PY_LONG_LONG, PyLong_AsLongLong(x)) +#endif } } { @@ -2763,8 +3128,10 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, #endif if (sizeof(int) <= sizeof(unsigned long)) { __PYX_VERIFY_RETURN_INT_EXC(int, unsigned long, PyLong_AsUnsignedLong(x)) +#ifdef HAVE_LONG_LONG } else if (sizeof(int) <= sizeof(unsigned PY_LONG_LONG)) { __PYX_VERIFY_RETURN_INT_EXC(int, unsigned PY_LONG_LONG, PyLong_AsUnsignedLongLong(x)) +#endif } } else { #if CYTHON_USE_PYLONG_INTERNALS @@ -2831,8 +3198,10 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, #endif if (sizeof(int) <= sizeof(long)) { __PYX_VERIFY_RETURN_INT_EXC(int, long, PyLong_AsLong(x)) +#ifdef HAVE_LONG_LONG } else if (sizeof(int) <= sizeof(PY_LONG_LONG)) { __PYX_VERIFY_RETURN_INT_EXC(int, PY_LONG_LONG, PyLong_AsLongLong(x)) +#endif } } { @@ -2948,8 +3317,10 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, #endif if (sizeof(long) <= sizeof(unsigned long)) { __PYX_VERIFY_RETURN_INT_EXC(long, unsigned long, PyLong_AsUnsignedLong(x)) +#ifdef HAVE_LONG_LONG } else if (sizeof(long) <= sizeof(unsigned PY_LONG_LONG)) { __PYX_VERIFY_RETURN_INT_EXC(long, unsigned PY_LONG_LONG, PyLong_AsUnsignedLongLong(x)) +#endif } } else { #if CYTHON_USE_PYLONG_INTERNALS @@ -3016,8 +3387,10 @@ static void __Pyx_AddTraceback(const char *funcname, int c_line, #endif if (sizeof(long) <= sizeof(long)) { __PYX_VERIFY_RETURN_INT_EXC(long, long, PyLong_AsLong(x)) +#ifdef HAVE_LONG_LONG } else if (sizeof(long) <= sizeof(PY_LONG_LONG)) { __PYX_VERIFY_RETURN_INT_EXC(long, PY_LONG_LONG, PyLong_AsLongLong(x)) +#endif } } { @@ -3183,7 +3556,9 @@ static CYTHON_INLINE int __Pyx_PyObject_IsTrue(PyObject* x) { else return PyObject_IsTrue(x); } static CYTHON_INLINE PyObject* __Pyx_PyNumber_IntOrLong(PyObject* x) { +#if CYTHON_USE_TYPE_SLOTS PyNumberMethods *m; +#endif const char *name = NULL; PyObject *res = NULL; #if PY_MAJOR_VERSION < 3 @@ -3192,8 +3567,9 @@ static CYTHON_INLINE PyObject* __Pyx_PyNumber_IntOrLong(PyObject* x) { if (PyLong_Check(x)) #endif return __Pyx_NewRef(x); +#if CYTHON_USE_TYPE_SLOTS m = Py_TYPE(x)->tp_as_number; -#if PY_MAJOR_VERSION < 3 + #if PY_MAJOR_VERSION < 3 if (m && m->nb_int) { name = "int"; res = PyNumber_Int(x); @@ -3202,11 +3578,14 @@ static CYTHON_INLINE PyObject* __Pyx_PyNumber_IntOrLong(PyObject* x) { name = "long"; res = PyNumber_Long(x); } -#else + #else if (m && m->nb_int) { name = "int"; res = PyNumber_Long(x); } + #endif +#else + res = PyNumber_Int(x); #endif if (res) { #if PY_MAJOR_VERSION < 3 diff --git a/tests/base/test_cyl.py b/tests/base/test_cyl.py index badea541b..810665069 100644 --- a/tests/base/test_cyl.py +++ b/tests/base/test_cyl.py @@ -204,7 +204,11 @@ def test_getInterpMatCartMesh_Faces(self): def test_getInterpMatCartMesh_Faces2Edges(self): Mr = discretize.TensorMesh([100, 100, 2], x0='CC0') - Mc = discretize.CylMesh([np.ones(10)/5, 1, 10], x0='0C0', cartesianOrigin=[-0.2, -0.2, 0]) + Mc = discretize.CylMesh( + [np.ones(10)/5, 1, 10], x0='0C0', cartesianOrigin=[-0.2, -0.2, 0] + ) + + self.assertTrue((Mc.cartesianOrigin == [-0.2, -0.2, 0]).all()) Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E') mf = np.ones(Mc.nF) From d9d59cc991bf4130cf45244125e5339198e505e7 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 09:38:19 -0700 Subject: [PATCH 11/34] pass x0 as a kwarg --- discretize/MeshIO.py | 6 +++--- tests/base/test_tensor.py | 37 +++++++++++++++++++++++++------------ 2 files changed, 28 insertions(+), 15 deletions(-) diff --git a/discretize/MeshIO.py b/discretize/MeshIO.py index 75a0ef7ea..31eec2022 100644 --- a/discretize/MeshIO.py +++ b/discretize/MeshIO.py @@ -91,7 +91,7 @@ def unpackdx(fid, nrows): z0 = z0 - sum(dz) dz = dz[::-1] # Make the mesh - tensMsh = TensorMesh([dx, dz], (x0, z0)) + tensMsh = TensorMesh([dx, dz], x0=(x0, z0)) fopen.close() @@ -161,7 +161,7 @@ def readVTK(TensorMesh, fileName): x0 = np.array([xR, yR, zR]) # Make the object - tensMsh = TensorMesh([hx, hy, hz], x0) + tensMsh = TensorMesh([hx, hy, hz], x0=x0) # Grap the models models = {} @@ -464,7 +464,7 @@ def readUBC(TreeMesh, meshFile): simpegPointers = np.concatenate((simpegCellPt, simpegLevel.reshape((-1, 1))), axis=1) # Make the tree mesh - mesh = TreeMesh([h1, h2, h3], x0) + mesh = TreeMesh([h1, h2, h3], x0=x0) mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) # Figure out the reordering diff --git a/tests/base/test_tensor.py b/tests/base/test_tensor.py index 1da60e7f1..3f1ab697d 100644 --- a/tests/base/test_tensor.py +++ b/tests/base/test_tensor.py @@ -13,7 +13,7 @@ def setUp(self): a = np.array([1, 1, 1]) b = np.array([1, 2]) c = np.array([1, 4]) - self.mesh2 = discretize.TensorMesh([a, b], [3, 5]) + self.mesh2 = discretize.TensorMesh([a, b], x0=[3, 5]) self.mesh3 = discretize.TensorMesh([a, b, c]) def test_vectorN_2D(self): @@ -33,12 +33,18 @@ def test_vectorCC_2D(self): self.assertTrue(xtest and ytest) def test_area_3D(self): - test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) + test_area = np.array([ + 1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, + 2, 2, 1, 1, 1, 2, 2, 2 + ]) t1 = np.all(self.mesh3.area == test_area) self.assertTrue(t1) def test_vol_3D(self): - test_vol = np.array([1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8]) + test_vol = np.array([ + 1, 1, 1, 2, 2, 2, 4, 4, 4, 8, 8, 8 + ]) t1 = np.all(self.mesh3.vol == test_vol) self.assertTrue(t1) @@ -48,12 +54,19 @@ def test_vol_2D(self): self.assertTrue(t1) def test_edge_3D(self): - test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]) + test_edge = np.array([ + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, + 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4 + ]) t1 = np.all(self.mesh3.edge == test_edge) self.assertTrue(t1) def test_edge_2D(self): - test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2]) + test_edge = np.array([ + 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2 + ]) t1 = np.all(self.mesh2.edge == test_edge) self.assertTrue(t1) @@ -68,24 +81,24 @@ def test_printing(self): print(discretize.TensorMesh([10, 10, 10])) def test_centering(self): - M1d = discretize.TensorMesh([10], 'C') - M2d = discretize.TensorMesh([10, 10], 'CC') - M3d = discretize.TensorMesh([10, 10, 10], 'CCC') + M1d = discretize.TensorMesh([10], x0='C') + M2d = discretize.TensorMesh([10, 10], x0='CC') + M3d = discretize.TensorMesh([10, 10, 10], x0='CCC') self.assertLess(np.abs(M1d.x0 + 0.5).sum(), TOL) self.assertLess(np.abs(M2d.x0 + 0.5).sum(), TOL) self.assertLess(np.abs(M3d.x0 + 0.5).sum(), TOL) def test_negative(self): - M1d = discretize.TensorMesh([10], 'N') + M1d = discretize.TensorMesh([10], x0='N') self.assertRaises(Exception, discretize.TensorMesh, [10], 'F') - M2d = discretize.TensorMesh([10, 10], 'NN') - M3d = discretize.TensorMesh([10, 10, 10], 'NNN') + M2d = discretize.TensorMesh([10, 10], x0='NN') + M3d = discretize.TensorMesh([10, 10, 10], x0='NNN') self.assertLess(np.abs(M1d.x0 + 1.0).sum(), TOL) self.assertLess(np.abs(M2d.x0 + 1.0).sum(), TOL) self.assertLess(np.abs(M3d.x0 + 1.0).sum(), TOL) def test_cent_neg(self): - M3d = discretize.TensorMesh([10, 10, 10], 'C0N') + M3d = discretize.TensorMesh([10, 10, 10], x0='C0N') self.assertLess(np.abs(M3d.x0 + np.r_[0.5, 0, 1.0]).sum(), TOL) def test_tensor(self): From 2606404d80b6180049263642d4dc06d06c617cb7 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 10:02:53 -0700 Subject: [PATCH 12/34] make levels a property in tree mesh. Pass x0 as a kwarg --- discretize/TreeMesh.py | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index 8c7be0397..92c52e04b 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -81,6 +81,7 @@ import numpy as np import scipy.sparse as sp +import properties from discretize import utils from . import TreeUtils @@ -98,24 +99,38 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): _meshType = 'TREE' - def __init__(self, h, x0=None, levels=None): + levels = properties.Integer( + "discretization level", + min=0 + ) + + def __init__(self, h, **kwargs): assert type(h) is list, 'h must be a list' assert len(h) in [2, 3], "TreeMesh is only in 2D or 3D." - BaseTensorMesh.__init__(self, h, x0) + if 'levels' in kwargs.keys(): + self.levels = kwargs.pop('levels') + + BaseTensorMesh.__init__(self, h, **kwargs) - if levels is None: - levels = int(np.log2(len(self.h[0]))) - assert np.all(len(_) == 2**levels for _ in self.h), "must make h and levels match" + if self.levels is None: + self.levels = int(np.log2(len(self.h[0]))) - self._levels = levels - self._levelBits = int(np.ceil(np.sqrt(levels)))+1 + # self._levels = levels + self._levelBits = int(np.ceil(np.sqrt(self.levels)))+1 self.__dirty__ = True #: The numbering is dirty! self._cells = set() self._cells.add(0) + @properties.validator('levels') + def check_levels(self, change): + assert np.all( + len(_) == 2**change['value'] for _ in self.h + ), "must make h and levels match" + + @property def __dirty__(self): return (self.__dirtyFaces__ or @@ -147,9 +162,9 @@ def __dirty__(self, val): for p in deleteThese: if hasattr(self, p): delattr(self, p) - @property - def levels(self): - return self._levels + # @property + # def levels(self): + # return self._levels @property def fill(self): From 24accb014e6e86dbb674ef57abebf23f785c2c13 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 10:24:34 -0700 Subject: [PATCH 13/34] pass x0 as a kwarg in view --- discretize/View.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/discretize/View.py b/discretize/View.py index eb49f341d..d97d25e1a 100644 --- a/discretize/View.py +++ b/discretize/View.py @@ -291,7 +291,7 @@ def doSlice(v): if 'Z' not in normal: h2d.append(self.hz) x2d.append(self.x0[2]) - tM = self.__class__(h2d, x2d) #: Temp Mesh + tM = self.__class__(h2d, x0=x2d) #: Temp Mesh v2d = doSlice(v) if ax is None: @@ -405,7 +405,7 @@ def _plotImage2D( hy = np.ones(nyi)*self.hy.sum()/nyi x0_y = self.x0[1] - tMi = self.__class__([hx, hy], np.r_[x0_x, x0_y]) + tMi = self.__class__([hx, hy], x0=np.r_[x0_x, x0_y]) P = self.getInterpolationMat(tMi.gridCC, 'CC', zerosOutside=True) Ui = tMi.r(P*mkvc(U), 'CC', 'CC', 'M') From d347514207c3cfe87b330c132d6f6f8127aa5088 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 10:25:03 -0700 Subject: [PATCH 14/34] add properties docs to intersphinx, remove quantified code from readme --- README.rst | 3 --- docs/conf.py | 3 ++- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/README.rst b/README.rst index 8dd96b8e9..d38a9a405 100644 --- a/README.rst +++ b/README.rst @@ -21,9 +21,6 @@ discretize :target: https://www.codacy.com/app/lindseyheagy/discretize?utm_source=github.com&utm_medium=referral&utm_content=simpeg/discretize&utm_campaign=Badge_Grade :alt: codacy status -.. image:: https://www.quantifiedcode.com/api/v1/project/97126041a7e7400fbf1c5555bed514dc/badge.svg - :target: https://www.quantifiedcode.com/app/project/97126041a7e7400fbf1c5555bed514dc - :alt: Code issues **discretize** - A python package for finite volume discretization. diff --git a/docs/conf.py b/docs/conf.py index 176de5ff1..35ea56686 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -254,7 +254,8 @@ 'python': ('https://docs.python.org/3/', None), 'numpy': ('https://docs.scipy.org/doc/numpy/', None), 'scipy': ('https://docs.scipy.org/doc/scipy/reference/', None), - 'matplotlib': ('http://matplotlib.org/', None) + 'matplotlib': ('http://matplotlib.org/', None), + 'properties': ('http://propertiespy.readthedocs.io/en/latest/', None), } From 7b03487ce8356abbcba3c27c91186d9f47765b80 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 13:18:07 -0700 Subject: [PATCH 15/34] remove requirement that x0 be passed as a kwarg --- discretize/BaseMesh.py | 4 ++-- discretize/CylMesh.py | 4 ++-- discretize/TensorMesh.py | 13 +++++++------ discretize/TreeMesh.py | 4 ++-- tests/base/test_cyl.py | 2 +- tests/base/test_tensor.py | 3 +-- 6 files changed, 15 insertions(+), 15 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index 000e3f782..deaddbe7a 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -27,11 +27,11 @@ class BaseMesh(properties.HasProperties): ) # Instantiate the class - def __init__(self, n, **kwargs): + def __init__(self, n, x0=None): self.n = n # number of dimensions # origin of the mesh - x0 = kwargs.pop('x0') + # x0 = kwargs.pop('x0') if x0 is None: self.x0 = np.zeros(len(self.n)) else: diff --git a/discretize/CylMesh.py b/discretize/CylMesh.py index da647e926..9129d9978 100644 --- a/discretize/CylMesh.py +++ b/discretize/CylMesh.py @@ -36,8 +36,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): shape=('*',) ) - def __init__(self, h_in=None, **kwargs): - BaseTensorMesh.__init__(self, h_in=h_in, **kwargs) + def __init__(self, h=None, x0=None, **kwargs): + BaseTensorMesh.__init__(self, h=h, x0=x0) assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" if self.dim == 2: print('Warning, a disk mesh has not been tested thoroughly.') diff --git a/discretize/TensorMesh.py b/discretize/TensorMesh.py index d0f62019e..f08787238 100644 --- a/discretize/TensorMesh.py +++ b/discretize/TensorMesh.py @@ -31,7 +31,10 @@ class BaseTensorMesh(BaseMesh): max_length=3 ) - def __init__(self, h_in=None, **kwargs): + def __init__(self, h=None, x0=None, **kwargs): + + h_in = h + x0_in = x0 # Cell widths if h_in is not None: @@ -68,9 +71,7 @@ def __init__(self, h_in=None, **kwargs): # Origin of the mesh x0 = np.zeros(len(h)) - if 'x0' in kwargs.keys(): - x0_in = kwargs.pop('x0') - + if x0_in is not None: assert len(h) == len(x0_in), "Dimension mismatch. x0 != len(h)" for i in range(len(h)): x_i, h_i = x0_in[i], h[i] @@ -543,8 +544,8 @@ class TensorMesh(BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, _meshType = 'TENSOR' - def __init__(self, h_in=None, **kwargs): - BaseTensorMesh.__init__(self, h_in, **kwargs) + def __init__(self, h=None, x0=None): + BaseTensorMesh.__init__(self, h=h, x0=x0) def __str__(self): outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index 92c52e04b..b235d40d6 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -104,14 +104,14 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): min=0 ) - def __init__(self, h, **kwargs): + def __init__(self, h, x0=None, **kwargs): assert type(h) is list, 'h must be a list' assert len(h) in [2, 3], "TreeMesh is only in 2D or 3D." if 'levels' in kwargs.keys(): self.levels = kwargs.pop('levels') - BaseTensorMesh.__init__(self, h, **kwargs) + BaseTensorMesh.__init__(self, h, x0, **kwargs) if self.levels is None: self.levels = int(np.log2(len(self.h[0]))) diff --git a/tests/base/test_cyl.py b/tests/base/test_cyl.py index 810665069..fc183a9b7 100644 --- a/tests/base/test_cyl.py +++ b/tests/base/test_cyl.py @@ -12,7 +12,7 @@ class TestCyl2DMesh(unittest.TestCase): def setUp(self): hx = np.r_[1, 1, 0.5] hz = np.r_[2, 1] - self.mesh = discretize.CylMesh([hx, 1, hz]) + self.mesh = discretize.CylMesh([hx, 1, hz], np.r_[0., 0., 0.]) def test_dim(self): self.assertTrue(self.mesh.dim == 3) diff --git a/tests/base/test_tensor.py b/tests/base/test_tensor.py index 3f1ab697d..2b01e0364 100644 --- a/tests/base/test_tensor.py +++ b/tests/base/test_tensor.py @@ -13,13 +13,12 @@ def setUp(self): a = np.array([1, 1, 1]) b = np.array([1, 2]) c = np.array([1, 4]) - self.mesh2 = discretize.TensorMesh([a, b], x0=[3, 5]) + self.mesh2 = discretize.TensorMesh([a, b], [3, 5]) self.mesh3 = discretize.TensorMesh([a, b, c]) def test_vectorN_2D(self): testNx = np.array([3, 4, 5, 6]) testNy = np.array([5, 6, 8]) - xtest = np.all(self.mesh2.vectorNx == testNx) ytest = np.all(self.mesh2.vectorNy == testNy) self.assertTrue(xtest and ytest) From f8c9c62a10aabfebdeb17ffb513875191d2366bd Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 17:46:56 -0700 Subject: [PATCH 16/34] check serialization of cyl mesh, allow kwargs in init of tensor meshes --- discretize/BaseMesh.py | 2 -- discretize/CurvilinearMesh.py | 4 ++- discretize/CylMesh.py | 2 +- discretize/TensorMesh.py | 15 ++++++--- discretize/TreeMesh.py | 52 ++++++++++++++++++----------- tests/base/test_properties.py | 61 ++++++++++++++++++++++++++++++++++- 6 files changed, 109 insertions(+), 27 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index deaddbe7a..c797b931c 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -30,8 +30,6 @@ class BaseMesh(properties.HasProperties): def __init__(self, n, x0=None): self.n = n # number of dimensions - # origin of the mesh - # x0 = kwargs.pop('x0') if x0 is None: self.x0 = np.zeros(len(self.n)) else: diff --git a/discretize/CurvilinearMesh.py b/discretize/CurvilinearMesh.py index a65f5b739..b3c0152e4 100644 --- a/discretize/CurvilinearMesh.py +++ b/discretize/CurvilinearMesh.py @@ -25,7 +25,9 @@ def normalize3D(x): return x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2)) -class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurviView): +class CurvilinearMesh( + BaseRectangularMesh, DiffOperators, InnerProducts, CurviView +): """CurvilinearMesh is a mesh class that deals with curvilinear meshes. Example of a curvilinear mesh: diff --git a/discretize/CylMesh.py b/discretize/CylMesh.py index 9129d9978..8c10f65c6 100644 --- a/discretize/CylMesh.py +++ b/discretize/CylMesh.py @@ -37,7 +37,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): ) def __init__(self, h=None, x0=None, **kwargs): - BaseTensorMesh.__init__(self, h=h, x0=x0) + BaseTensorMesh.__init__(self, h=h, x0=x0, **kwargs) assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" if self.dim == 2: print('Warning, a disk mesh has not been tested thoroughly.') diff --git a/discretize/TensorMesh.py b/discretize/TensorMesh.py index f08787238..0d615da53 100644 --- a/discretize/TensorMesh.py +++ b/discretize/TensorMesh.py @@ -89,8 +89,16 @@ def __init__(self, h=None, x0=None, **kwargs): "'C' to center, or 'N' to be negative.".format(i) ) + if 'n' in kwargs.keys(): + n = kwargs.pop('n') + assert (n == np.array([x.size for x in h])).all(), ( + "Dimension mismatch. The provided n doesn't " + ) + + n = np.array([x.size for x in h]) + super(BaseTensorMesh, self).__init__( - np.array([x.size for x in h]), x0=x0 + n, x0=x0 ) # Ensure h contains 1D vectors @@ -344,7 +352,6 @@ def getInterpolationMat(self, loc, locType='CC', zerosOutside=False): """ return self._getInterpolationMat(loc, locType, zerosOutside) - def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): """ Fast version of getFaceInnerProduct. @@ -544,8 +551,8 @@ class TensorMesh(BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, _meshType = 'TENSOR' - def __init__(self, h=None, x0=None): - BaseTensorMesh.__init__(self, h=h, x0=x0) + def __init__(self, h=None, x0=None, **kwargs): + BaseTensorMesh.__init__(self, h=h, x0=x0, **kwargs) def __str__(self): outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index b235d40d6..95cf51432 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -151,21 +151,19 @@ def __dirty__(self, val): self.__dirtySets__ = True deleteThese = [ - '__sortedCells', - '_gridCC', '_gridN', '_gridFx', '_gridFy', '_gridFz', '_gridEx', '_gridEy', '_gridEz', - '_area', '_edge', '_vol', - '_faceDiv', '_edgeCurl', '_nodalGrad', - '_aveFx2CC', '_aveFy2CC', '_aveFz2CC', '_aveF2CC', '_aveF2CCV', - '_aveEx2CC', '_aveEy2CC', '_aveEz2CC', '_aveE2CC', '_aveE2CCV', - '_aveN2CC', - ] + '__sortedCells', + '_gridCC', + '_gridN', '_gridFx', '_gridFy', '_gridFz', + '_gridEx', '_gridEy', '_gridEz', + '_area', '_edge', '_vol', + '_faceDiv', '_edgeCurl', '_nodalGrad', + '_aveFx2CC', '_aveFy2CC', '_aveFz2CC', '_aveF2CC', '_aveF2CCV', + '_aveEx2CC', '_aveEy2CC', '_aveEz2CC', '_aveE2CC', '_aveE2CCV', + '_aveN2CC', + ] for p in deleteThese: if hasattr(self, p): delattr(self, p) - # @property - # def levels(self): - # return self._levels - @property def fill(self): """How filled is the mesh compared to a TensorMesh? As a fraction: [0,1].""" @@ -432,7 +430,10 @@ def _pointer(self, index): def __contains__(self, v): return self._asIndex(v) in self._cells - def refine(self, function=None, recursive=True, cells=None, balance=True, verbose=False, _inRecursion=False): + def refine( + self, function=None, recursive=True, cells=None, balance=True, + verbose=False, _inRecursion=False + ): if type(function) in integer_types: level = function @@ -454,20 +455,29 @@ def refine(self, function=None, recursive=True, cells=None, balance=True, verbos elif type(result) in integer_types: do = result > p[-1] else: - raise Exception('You must tell the program what to refine. Use BOOL or INT (level)') + raise Exception( + 'You must tell the program what to refine. ' + 'Use BOOL or INT (level)' + ) if do: recurse += self._refineCell(cell, p) if verbose: print(' ', time.time() - tic) if recursive and len(recurse) > 0: - recurse += self.refine(function=function, recursive=True, cells=recurse, balance=balance, verbose=verbose, _inRecursion=True) + recurse += self.refine( + function=function, recursive=True, cells=recurse, + balance=balance, verbose=verbose, _inRecursion=True + ) if balance and not _inRecursion: self.balance() return recurse - def corsen(self, function=None, recursive=True, cells=None, balance=True, verbose=False, _inRecursion=False): + def corsen( + self, function=None, recursive=True, cells=None, balance=True, + verbose=False, _inRecursion=False + ): if type(function) in integer_types: level = function @@ -490,14 +500,20 @@ def corsen(self, function=None, recursive=True, cells=None, balance=True, verbos elif type(result) in integer_types: do = result < p[-1] else: - raise Exception('You must tell the program what to corsen. Use BOOL or INT (level)') + raise Exception( + 'You must tell the program what to corsen. Use BOOL or ' + 'INT (level)' + ) if do: recurse += self._corsenCell(cell, p) if verbose: print(' ', time.time() - tic) if recursive and len(recurse) > 0: - recurse += self.corsen(function=function, recursive=True, cells=recurse, balance=balance, verbose=verbose, _inRecursion=True) + recurse += self.corsen( + function=function, recursive=True, cells=recurse, + balance=balance, verbose=verbose, _inRecursion=True + ) if balance and not _inRecursion: self.balance() diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index 0c2e21c9d..baaaa7b80 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -17,7 +17,6 @@ def compare_meshes(mesh0, mesh1): 'mesh h[{}] different'.format(i) ) - class TensorTest(unittest.TestCase): n = [4, 5, 9] @@ -43,6 +42,66 @@ def test_copy(self): print('ok\n') +class CylTest(unittest.TestCase): + + n = [4, 1, 9] + + def setUp(self): + self.mesh = discretize.CylMesh(self.n, x0='00C') + + def test_save_load(self): + print('\nTesting save / load of Cyl Mesh ...') + mesh0 = self.mesh + f = mesh0.save() + mesh1 = discretize.utils.load_mesh(f) + compare_meshes(mesh0, mesh1) + os.remove(f) + print('ok\n') + + def test_copy(self): + print('\nTesting copy of Cyl Mesh ...') + mesh0 = self.mesh + mesh1 = mesh0.copy() + compare_meshes(mesh0, mesh1) + print('ok\n') + + +# class TreeTest(unittest.TestCase): + +# n = [4, 1, 9] + +# def setUp(self): +# M = discretize.TreeMesh([8, 8]) + +# def refine(cell): +# xyz = cell.center +# dist = ((xyz - [0.25, 0.25])**2).sum()**0.5 +# if dist < 0.25: +# return 3 +# return 2 + +# M.refine(refine) +# M.number() + +# self.mesh = M + +# def test_save_load(self): +# print('\nTesting save / load of Tree Mesh ...') +# mesh0 = self.mesh +# f = mesh0.save() +# mesh1 = discretize.utils.load_mesh(f) +# compare_meshes(mesh0, mesh1) +# os.remove(f) +# print('ok\n') + +# def test_copy(self): +# print('\nTesting copy of Tree Mesh ...') +# mesh0 = self.mesh +# mesh1 = mesh0.copy() +# compare_meshes(mesh0, mesh1) +# print('ok\n') + + if __name__ == '__main__': From 393851ff6b5abeb2b8e047e1f68afdbb6a0037e0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 18:04:39 -0700 Subject: [PATCH 17/34] make _cells a public variable so that it can be a property (see aranzgeo/properties#191, test serialization and deserialization of tree mesh --- discretize/TreeMesh.py | 51 +++++++++++++++++------------ tests/base/test_properties.py | 60 +++++++++++++++++------------------ 2 files changed, 61 insertions(+), 50 deletions(-) diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index 95cf51432..c7cbf8370 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -104,6 +104,15 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): min=0 ) + cells = properties.Set( + "cells in the tree mesh", + properties.Integer( + "cell number", + min=0 + ), + default=set() + ) + def __init__(self, h, x0=None, **kwargs): assert type(h) is list, 'h must be a list' assert len(h) in [2, 3], "TreeMesh is only in 2D or 3D." @@ -121,8 +130,10 @@ def __init__(self, h, x0=None, **kwargs): self.__dirty__ = True #: The numbering is dirty! - self._cells = set() - self._cells.add(0) + if 'cells' in kwargs.keys(): + self.cells = kwargs.pop('cells') + else: + self.cells.add(0) @properties.validator('levels') def check_levels(self, change): @@ -173,7 +184,7 @@ def fill(self): def maxLevel(self): """The maximum level used, which may be less than `levels`.""" l = 0 - for cell in self._cells: + for cell in self.cells: p = self._pointer(cell) l = max(l,p[-1]) return l @@ -221,7 +232,7 @@ def printH(hx, outStr=''): @property def nC(self): - return len(self._cells) + return len(self.cells) @property def nN(self): @@ -387,7 +398,7 @@ def ntEz(self): @property def _sortedCells(self): if getattr(self, '__sortedCells', None) is None: - self.__sortedCells = sorted(self._cells) + self.__sortedCells = sorted(self.cells) return self.__sortedCells @property @@ -428,7 +439,7 @@ def _pointer(self, index): return TreeUtils.point(self.dim, MAX_BITS, self._levelBits, index) def __contains__(self, v): - return self._asIndex(v) in self._cells + return self._asIndex(v) in self.cells def refine( self, function=None, recursive=True, cells=None, balance=True, @@ -443,7 +454,7 @@ def refine( self.__dirty__ = True if verbose: print('Refining Mesh') - cells = cells if cells is not None else sorted(self._cells) + cells = cells if cells is not None else sorted(self.cells) recurse = [] tic = time.time() for cell in cells: @@ -487,11 +498,11 @@ def corsen( self.__dirty__ = True if verbose: print('Corsening Mesh') - cells = cells if cells is not None else sorted(self._cells) + cells = cells if cells is not None else sorted(self.cells) recurse = [] tic = time.time() for cell in cells: - if cell not in self._cells: continue # already removed + if cell not in self.cells: continue # already removed p = self._pointer(cell) if p[-1] >= self.levels: continue result = function(Cell(self, cell, p)) @@ -526,8 +537,8 @@ def _refineCell(self, ind, pointer=None): raise IndexError(ind) children = self._childPointers(pointer, returnAll=True) for child in children: - self._cells.add(self._asIndex(child)) - self._cells.remove(ind) + self.cells.add(self._asIndex(child)) + self.cells.remove(ind) return [self._asIndex(child) for child in children] def _corsenCell(self, ind, pointer=None): @@ -538,9 +549,9 @@ def _corsenCell(self, ind, pointer=None): parent = self._parentPointer(pointer) children = self._childPointers(parent, returnAll=True) for child in children: - self._cells.remove(self._asIndex(child)) + self.cells.remove(self._asIndex(child)) parentInd = self._asIndex(parent) - self._cells.add(parentInd) + self.cells.add(parentInd) return [parentInd] def _asPointer(self, ind): @@ -665,7 +676,7 @@ def balance(self, recursive=True, cells=None, verbose=False, _inRecursion=False) self.__dirty__ = True if verbose: print('Balancing Mesh:') - cells = cells if cells is not None else sorted(self._cells) + cells = cells if cells is not None else sorted(self.cells) # calcDepth = lambda i: lambda A: i if type(A) is not list else max(map(calcDepth(i+1), A)) # flatten = lambda A: A if calcDepth(0)(A) == 1 else flatten([_ for __ in A for _ in (__ if type(__) is list else [__])]) @@ -705,7 +716,7 @@ def balance(self, recursive=True, cells=None, verbose=False, _inRecursion=False) @property def gridCC(self): if getattr(self, '_gridCC', None) is None: - self._gridCC = np.zeros((len(self._cells),self.dim)) + self._gridCC = np.zeros((len(self.cells),self.dim)) for ii, ind in enumerate(self._sortedCells): p = self._asPointer(ind) self._gridCC[ii, :] = self._cellC(p) + self.x0 @@ -764,7 +775,7 @@ def gridEz(self): @property def vol(self): if getattr(self, '_vol', None) is None: - self._vol = np.zeros(len(self._cells)) + self._vol = np.zeros(len(self.cells)) for ii, ind in enumerate(self._sortedCells): p = self._asPointer(ind) self._vol[ii] = np.prod(self._cellH(p)) @@ -807,7 +818,7 @@ def _createNumberingSets(self, force=False): self._edgesZ = set() - for ind in self._cells: + for ind in self.cells: p = self._asPointer(ind) w = self._levelWidth(p[-1]) if self.dim == 2: @@ -878,7 +889,7 @@ def _numberCells(self, force=False): return self._cc2i = dict() self._i2cc = dict() - for ii, c in enumerate(sorted(self._cells)): + for ii, c in enumerate(sorted(self.cells)): self._cc2i[c] = ii self._i2cc[ii] = c self.__dirtyCells__ = False @@ -901,7 +912,7 @@ def _numberFaces(self, force=False): return self._createNumberingSets(force=force) - for ind in self._cells: + for ind in self.cells: p = self._asPointer(ind) w = self._levelWidth(p[-1]) @@ -2273,7 +2284,7 @@ def doSlice(v): _ = antiNormalInd X = [] Y = [] - for cell in self._cells: + for cell in self.cells: p = self._pointer(cell) n, h = self._cellN(p), self._cellH(p) if n[normalInd]indLoc: diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index baaaa7b80..f1fd02d11 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -17,6 +17,10 @@ def compare_meshes(mesh0, mesh1): 'mesh h[{}] different'.format(i) ) + if hasattr(mesh0, 'cells'): + assert (mesh0.cells == mesh1.cells) + + class TensorTest(unittest.TestCase): n = [4, 5, 9] @@ -66,42 +70,38 @@ def test_copy(self): print('ok\n') -# class TreeTest(unittest.TestCase): - -# n = [4, 1, 9] - -# def setUp(self): -# M = discretize.TreeMesh([8, 8]) - -# def refine(cell): -# xyz = cell.center -# dist = ((xyz - [0.25, 0.25])**2).sum()**0.5 -# if dist < 0.25: -# return 3 -# return 2 +class TreeTest(unittest.TestCase): -# M.refine(refine) -# M.number() + def setUp(self): + M = discretize.TreeMesh([8, 8]) -# self.mesh = M + def refine(cell): + xyz = cell.center + dist = ((xyz - [0.25, 0.25])**2).sum()**0.5 + if dist < 0.25: + return 3 + return 2 -# def test_save_load(self): -# print('\nTesting save / load of Tree Mesh ...') -# mesh0 = self.mesh -# f = mesh0.save() -# mesh1 = discretize.utils.load_mesh(f) -# compare_meshes(mesh0, mesh1) -# os.remove(f) -# print('ok\n') + M.refine(refine) + M.number() -# def test_copy(self): -# print('\nTesting copy of Tree Mesh ...') -# mesh0 = self.mesh -# mesh1 = mesh0.copy() -# compare_meshes(mesh0, mesh1) -# print('ok\n') + self.mesh = M + def test_save_load(self): + print('\nTesting save / load of Tree Mesh ...') + mesh0 = self.mesh + f = mesh0.save() + mesh1 = discretize.utils.load_mesh(f) + compare_meshes(mesh0, mesh1) + os.remove(f) + print('ok\n') + def test_copy(self): + print('\nTesting copy of Tree Mesh ...') + mesh0 = self.mesh + mesh1 = mesh0.copy() + compare_meshes(mesh0, mesh1) + print('ok\n') if __name__ == '__main__': From e770ff9f7f8c522e7842475f5f1c9257597821f0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 1 Sep 2017 18:12:56 -0700 Subject: [PATCH 18/34] add tests to check edges, volumes and areas are still the same --- discretize/TreeMesh.py | 2 +- tests/base/test_properties.py | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index c7cbf8370..be6b66580 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -186,7 +186,7 @@ def maxLevel(self): l = 0 for cell in self.cells: p = self._pointer(cell) - l = max(l,p[-1]) + l = max(l, p[-1]) return l def __str__(self): diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index f1fd02d11..e6ab3c524 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -1,9 +1,12 @@ import unittest import os +import numpy as np import discretize def compare_meshes(mesh0, mesh1): + + # check some basic properties assert mesh0.nC == mesh1.nC , ( 'Number of cells not the same, {} != {}'.format(mesh0.nC, mesh1.nC) ) @@ -17,6 +20,11 @@ def compare_meshes(mesh0, mesh1): 'mesh h[{}] different'.format(i) ) + # check edges, faces, volumes + assert (mesh0.edge == mesh1.edge).all() + assert (mesh0.area == mesh1.area).all() + assert (mesh0.vol == mesh1.vol).all() + if hasattr(mesh0, 'cells'): assert (mesh0.cells == mesh1.cells) From 8b061fecbda0f5cc76250ca0539a0c632fb62661 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 2 Sep 2017 10:18:06 -0700 Subject: [PATCH 19/34] add a subclass of the properties Array to overcome the size issue discussed in aranzgeo/properties#192 --- discretize/CurvilinearMesh.py | 105 +++++++++++++++++++++++++++++----- tests/base/test_properties.py | 46 +++++++++++++-- 2 files changed, 132 insertions(+), 19 deletions(-) diff --git a/discretize/CurvilinearMesh.py b/discretize/CurvilinearMesh.py index b3c0152e4..fa8eb71ce 100644 --- a/discretize/CurvilinearMesh.py +++ b/discretize/CurvilinearMesh.py @@ -1,5 +1,7 @@ from __future__ import print_function import numpy as np +import properties +from properties.math import TYPE_MAPPINGS from discretize import utils from discretize.BaseMesh import BaseRectangularMesh @@ -25,6 +27,59 @@ def normalize3D(x): return x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2)) +class Array(properties.Array): + + @property + def shape(self): + """Valid array shape. + + Must be a tuple with integer or '*' engries corresponding to valid + array shapes. '*' means the dimension can be any length. + """ + return getattr(self, '_shape', '*') + + @shape.setter + def shape(self, value): + if not isinstance(value, tuple): + if not (isinstance(value, str) and value == '*'): + raise TypeError( + "{}: Invalid shape - must be a tuple (e.g. ('*',3) for an " + "array of length-3 arrays) or '*'".format(value) + ) + else: + for shp in value: + if shp != '*' and not isinstance(shp, integer_types): + raise TypeError( + "{}: Invalid shape - values must be '*' or int".format( + value + ) + ) + self._shape = value + + def validate(self, instance, value): + """Determine if array is valid based on shape and dtype""" + if not isinstance(value, (tuple, list, np.ndarray)): + self.error(instance, value) + value = self.wrapper(value) + if not isinstance(value, np.ndarray): + raise NotImplementedError( + 'Array validation is only implmented for wrappers that are ' + 'subclasses of numpy.ndarray' + ) + if value.dtype.kind not in (TYPE_MAPPINGS[typ] for typ in self.dtype): + self.error(instance, value) + if not isinstance(self.shape, str): + if len(self.shape) != value.ndim: + self.error(instance, value) + for i, shp in enumerate(self.shape): + if shp != '*' and value.shape[i] != shp: + self.error(instance, value) + else: + if self.shape != '*': + self.error(isinstance, value) + return value + + class CurvilinearMesh( BaseRectangularMesh, DiffOperators, InnerProducts, CurviView ): @@ -43,28 +98,48 @@ class CurvilinearMesh( _meshType = 'Curv' - def __init__(self, nodes): - assert type(nodes) == list, ("'nodes' variable must be a list of " - "np.ndarray") - assert len(nodes) > 1, "len(node) must be greater than 1" + nodes = properties.List( + "List of arrays describing the node locations", + Array( + "node locations in a dimension", + shape='*' + ), + ) - for i, nodes_i in enumerate(nodes): - assert isinstance(nodes_i, np.ndarray), ("nodes[{0:d}] is not a" - "numpy array.".format(i)) - assert nodes_i.shape == nodes[0].shape, ("nodes[{0:d}] is not the " - "same shape as nodes[0]" - .format(i)) + def __init__(self, nodes=nodes, **kwargs): - assert len(nodes[0].shape) == len(nodes), "Dimension mismatch" - assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh." + self.nodes = nodes - BaseRectangularMesh.__init__(self, np.array(nodes[0].shape)-1, None) + if 'n' in kwargs.keys(): + n = kwargs.pop('n') + assert (n == np.array(self.nodes[0].shape)-1).all(), ( + "Unexpected n-values. {} was provided, {} was expected".format( + n, np.array(self.nodes[0].shape)-1 + ) + ) + else: + n = np.array(self.nodes[0].shape)-1 + + BaseRectangularMesh.__init__(self, n, **kwargs) # Save nodes to private variable _gridN as vectors - self._gridN = np.ones((nodes[0].size, self.dim)) - for i, node_i in enumerate(nodes): + self._gridN = np.ones((self.nodes[0].size, self.dim)) + for i, node_i in enumerate(self.nodes): self._gridN[:, i] = utils.mkvc(node_i.astype(float)) + @properties.validator('nodes') + def check_nodes(self, change): + assert len(change['value']) > 1, "len(node) must be greater than 1" + + for i, change['value'][i] in enumerate(change['value']): + assert change['value'][i].shape == change['value'][0].shape, ( + "change['value'][{0:d}] is not the same shape as change['value'][0]".format(i) + ) + + assert len(change['value'][0].shape) == len(change['value']), "Dimension mismatch" + assert len(change['value'][0].shape) > 1, "Not worth using Curv for a 1D mesh." + + @property def gridCC(self): """ diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index e6ab3c524..1c3fcabb0 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -15,19 +15,26 @@ def compare_meshes(mesh0, mesh1): 'x0 different. {} != {}'.format(mesh0.x0, mesh1.x0) ) - for i in range(len(mesh0.h)): - assert (mesh0.h[i] == mesh1.h[i]).all(), ( - 'mesh h[{}] different'.format(i) - ) + if hasattr(mesh0, 'h'): + for i in range(len(mesh0.h)): + assert (mesh0.h[i] == mesh1.h[i]).all(), ( + 'mesh h[{}] different'.format(i) + ) # check edges, faces, volumes assert (mesh0.edge == mesh1.edge).all() assert (mesh0.area == mesh1.area).all() assert (mesh0.vol == mesh1.vol).all() + # Tree mesh specific if hasattr(mesh0, 'cells'): assert (mesh0.cells == mesh1.cells) + # curvi-specific + if hasattr(mesh0, 'nodes'): + for i in range(len(mesh0.nodes)): + assert (mesh0.nodes[i] == mesh1.nodes[i]).all() + class TensorTest(unittest.TestCase): @@ -112,5 +119,36 @@ def test_copy(self): print('ok\n') +class CurviTest(unittest.TestCase): + + def setUp(self): + a = np.array([1, 1, 1]) + b = np.array([1, 2]) + c = np.array([1, 4]) + + def gridIt(h): return [np.cumsum(np.r_[0, x]) for x in h] + + X, Y, Z = discretize.utils.ndgrid(gridIt([a, b, c]), vector=False) + self.mesh = discretize.CurvilinearMesh([X, Y, Z]) + + def test_save_load(self): + print('\nTesting save / load of Curvi Mesh ...') + mesh0 = self.mesh + f = mesh0.save() + mesh1 = discretize.utils.load_mesh(f) + compare_meshes(mesh0, mesh1) + os.remove(f) + print('ok\n') + + def test_copy(self): + print('\nTesting copy of Curvi Mesh ...') + mesh0 = self.mesh + mesh1 = mesh0.copy() + compare_meshes(mesh0, mesh1) + print('ok\n') + + + + if __name__ == '__main__': unittest.main() From f471d17f54d14298a236d22977dbe40c6cfdc708 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 2 Sep 2017 10:53:23 -0700 Subject: [PATCH 20/34] add docstring for Array. add nitpick ignore for docs --- discretize/CurvilinearMesh.py | 17 +++++++++++++++++ docs/conf.py | 4 ++++ 2 files changed, 21 insertions(+) diff --git a/discretize/CurvilinearMesh.py b/discretize/CurvilinearMesh.py index fa8eb71ce..04b6d1c23 100644 --- a/discretize/CurvilinearMesh.py +++ b/discretize/CurvilinearMesh.py @@ -28,6 +28,23 @@ def normalize3D(x): class Array(properties.Array): + """ + An Array property (inheriting from + :class:`properties Array ` that allows for arbitrary + shape. + + **Available keywords** (in addition to those inherited from + :ref:`Property `): + + * **shape** - '*' or Tuple that describes the allowed shape of the array. + Length of shape tuple corresponds to number of dimensions; values + correspond to the allowed length for each dimension. These values + may be integers or '*' for any length. + For example, an n x 3 array would be shape ('*', 3). + The default value is ('*',). '*' indicated an arbitrary shape. + * **dtype** - Allowed data type for the array. May be float, int, + bool, or a tuple containing any of these. The default is (float, int). + """ @property def shape(self): diff --git a/docs/conf.py b/docs/conf.py index 35ea56686..9459bc10e 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -302,3 +302,7 @@ def _supress_nonlocal_image_warn(self, msg, node, **kwargs): self._warnfunc(msg, '{0!s}:{1!s}'.format(*get_source_line(node))) supress_nonlocal_image_warn() + +nitpick_ignore = [ + ('py:class', 'discretize.CurvilinearMesh.Array') +] From fd85204c0b9afc147b48ec8e5eda245ca37cac17 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 2 Sep 2017 11:30:15 -0700 Subject: [PATCH 21/34] set in the meshIO rather than for reading in UBC meshes --- discretize/MeshIO.py | 2 +- tests/base/test_properties.py | 6 ++++++ tests/tree/test_tree_io.py | 2 +- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/discretize/MeshIO.py b/discretize/MeshIO.py index 31eec2022..25ce4717a 100644 --- a/discretize/MeshIO.py +++ b/discretize/MeshIO.py @@ -465,7 +465,7 @@ def readUBC(TreeMesh, meshFile): # Make the tree mesh mesh = TreeMesh([h1, h2, h3], x0=x0) - mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) + mesh.cells = set([mesh._index(p) for p in simpegPointers.tolist()]) # Figure out the reordering mesh._simpegReorderUBC = np.argsort( diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index 1c3fcabb0..008098fe6 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -26,6 +26,12 @@ def compare_meshes(mesh0, mesh1): assert (mesh0.area == mesh1.area).all() assert (mesh0.vol == mesh1.vol).all() + # check strings + # print(mesh0.__str__()) + # print(mesh1.__str__()) + print("{}".format(mesh0.__str__())) + print("{}".format(mesh1.__str__())) + # Tree mesh specific if hasattr(mesh0, 'cells'): assert (mesh0.cells == mesh1.cells) diff --git a/tests/tree/test_tree_io.py b/tests/tree/test_tree_io.py index 08f81e9e3..2d7b321cd 100644 --- a/tests/tree/test_tree_io.py +++ b/tests/tree/test_tree_io.py @@ -32,7 +32,7 @@ def test_UBCfiles(self): meshUBC = discretize.TreeMesh.readUBC('temp.msh') vecUBC = meshUBC.readModelUBC('arange.txt') - # The mesh + assert(mesh.nC == meshUBC.nC) assert mesh.__str__() == meshUBC.__str__() assert np.sum(mesh.gridCC - meshUBC.gridCC) == 0 assert np.sum(vec - vecUBC) == 0 From abee61908063ab842100f9dec45e4436d7d74713 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 2 Sep 2017 11:41:51 -0700 Subject: [PATCH 22/34] add a few more property checks in the testing --- tests/base/test_properties.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index 008098fe6..6b795039e 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -15,6 +15,19 @@ def compare_meshes(mesh0, mesh1): 'x0 different. {} != {}'.format(mesh0.x0, mesh1.x0) ) + assert (mesh0.nE == mesh1.nE) + assert (mesh0.nF == mesh1.nF) + assert (mesh0.nN == mesh1.nN) + + if hasattr(mesh0, 'vnC'): + assert (mesh0.vnC == mesh1.vnC).all() + assert (mesh0.vnE == mesh1.vnE).all() + assert (mesh0.vnF == mesh1.vnF).all() + assert (mesh0.vnN == mesh1.vnN).all() + + assert (mesh0.normals == mesh1.normals).all() + assert (mesh0.tangents == mesh1.tangents).all() + if hasattr(mesh0, 'h'): for i in range(len(mesh0.h)): assert (mesh0.h[i] == mesh1.h[i]).all(), ( @@ -26,12 +39,6 @@ def compare_meshes(mesh0, mesh1): assert (mesh0.area == mesh1.area).all() assert (mesh0.vol == mesh1.vol).all() - # check strings - # print(mesh0.__str__()) - # print(mesh1.__str__()) - print("{}".format(mesh0.__str__())) - print("{}".format(mesh1.__str__())) - # Tree mesh specific if hasattr(mesh0, 'cells'): assert (mesh0.cells == mesh1.cells) From 5507fc101f2f4fb5ed9cc3666a9b2b295aedccf2 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 5 Sep 2017 12:46:00 -0700 Subject: [PATCH 23/34] use the union property for an arbitrary array size --- discretize/CurvilinearMesh.py | 98 ++++++----------------------------- 1 file changed, 17 insertions(+), 81 deletions(-) diff --git a/discretize/CurvilinearMesh.py b/discretize/CurvilinearMesh.py index 04b6d1c23..be1d7d5e4 100644 --- a/discretize/CurvilinearMesh.py +++ b/discretize/CurvilinearMesh.py @@ -27,76 +27,6 @@ def normalize3D(x): return x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2)) -class Array(properties.Array): - """ - An Array property (inheriting from - :class:`properties Array ` that allows for arbitrary - shape. - - **Available keywords** (in addition to those inherited from - :ref:`Property `): - - * **shape** - '*' or Tuple that describes the allowed shape of the array. - Length of shape tuple corresponds to number of dimensions; values - correspond to the allowed length for each dimension. These values - may be integers or '*' for any length. - For example, an n x 3 array would be shape ('*', 3). - The default value is ('*',). '*' indicated an arbitrary shape. - * **dtype** - Allowed data type for the array. May be float, int, - bool, or a tuple containing any of these. The default is (float, int). - """ - - @property - def shape(self): - """Valid array shape. - - Must be a tuple with integer or '*' engries corresponding to valid - array shapes. '*' means the dimension can be any length. - """ - return getattr(self, '_shape', '*') - - @shape.setter - def shape(self, value): - if not isinstance(value, tuple): - if not (isinstance(value, str) and value == '*'): - raise TypeError( - "{}: Invalid shape - must be a tuple (e.g. ('*',3) for an " - "array of length-3 arrays) or '*'".format(value) - ) - else: - for shp in value: - if shp != '*' and not isinstance(shp, integer_types): - raise TypeError( - "{}: Invalid shape - values must be '*' or int".format( - value - ) - ) - self._shape = value - - def validate(self, instance, value): - """Determine if array is valid based on shape and dtype""" - if not isinstance(value, (tuple, list, np.ndarray)): - self.error(instance, value) - value = self.wrapper(value) - if not isinstance(value, np.ndarray): - raise NotImplementedError( - 'Array validation is only implmented for wrappers that are ' - 'subclasses of numpy.ndarray' - ) - if value.dtype.kind not in (TYPE_MAPPINGS[typ] for typ in self.dtype): - self.error(instance, value) - if not isinstance(self.shape, str): - if len(self.shape) != value.ndim: - self.error(instance, value) - for i, shp in enumerate(self.shape): - if shp != '*' and value.shape[i] != shp: - self.error(instance, value) - else: - if self.shape != '*': - self.error(isinstance, value) - return value - - class CurvilinearMesh( BaseRectangularMesh, DiffOperators, InnerProducts, CurviView ): @@ -117,13 +47,18 @@ class CurvilinearMesh( nodes = properties.List( "List of arrays describing the node locations", - Array( - "node locations in a dimension", - shape='*' + prop=properties.Union( + "node locations in an n-dimensional array", + props=[ + properties.Array('2D array', shape=('*', '*')), + properties.Array('3D array', shape=('*', '*', '*')), + ] ), + min_length=2, + max_length=3 ) - def __init__(self, nodes=nodes, **kwargs): + def __init__(self, nodes=None, **kwargs): self.nodes = nodes @@ -150,12 +85,13 @@ def check_nodes(self, change): for i, change['value'][i] in enumerate(change['value']): assert change['value'][i].shape == change['value'][0].shape, ( - "change['value'][{0:d}] is not the same shape as change['value'][0]".format(i) + "change['value'][{0:d}] is not the same shape as " + "change['value'][0]".format(i) ) - assert len(change['value'][0].shape) == len(change['value']), "Dimension mismatch" - assert len(change['value'][0].shape) > 1, "Not worth using Curv for a 1D mesh." - + assert len(change['value'][0].shape) == len(change['value']), ( + "Dimension mismatch" + ) @property def gridCC(self): @@ -163,9 +99,9 @@ def gridCC(self): Cell-centered grid """ if getattr(self, '_gridCC', None) is None: - self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:, i] - for i in range(self.dim)]).reshape( - (-1, self.dim), order='F') + self._gridCC = np.concatenate( + [self.aveN2CC*self.gridN[:, i] for i in range(self.dim)] + ).reshape((-1, self.dim), order='F') return self._gridCC @property From 29305036c7c613fc8e34f5f563dc7f2ed97fc99a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 5 Sep 2017 12:54:04 -0700 Subject: [PATCH 24/34] whitespace cleanup in TensorMesh --- discretize/TensorMesh.py | 215 ++++++++++++++++++++++++--------------- 1 file changed, 135 insertions(+), 80 deletions(-) diff --git a/discretize/TensorMesh.py b/discretize/TensorMesh.py index 0d615da53..2cf856cd8 100644 --- a/discretize/TensorMesh.py +++ b/discretize/TensorMesh.py @@ -127,12 +127,16 @@ def vectorNx(self): @property def vectorNy(self): """Nodal grid vector (1D) in the y direction.""" - return None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1] + return ( + None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1] + ) @property def vectorNz(self): """Nodal grid vector (1D) in the z direction.""" - return None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2] + return ( + None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2] + ) @property def vectorCCx(self): @@ -142,12 +146,18 @@ def vectorCCx(self): @property def vectorCCy(self): """Cell-centered grid vector (1D) in the y direction.""" - return None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1] + return ( + None if self.dim < 2 else + np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1] + ) @property def vectorCCz(self): """Cell-centered grid vector (1D) in the z direction.""" - return None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2] + return ( + None if self.dim < 3 else + np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2] + ) @property def gridCC(self): @@ -162,37 +172,43 @@ def gridN(self): @property def gridFx(self): """Face staggered grid in the x direction.""" - if self.nFx == 0: return + if self.nFx == 0: + return return self._getTensorGrid('Fx') @property def gridFy(self): """Face staggered grid in the y direction.""" - if self.nFy == 0 or self.dim < 2: return + if self.nFy == 0 or self.dim < 2: + return return self._getTensorGrid('Fy') @property def gridFz(self): """Face staggered grid in the z direction.""" - if self.nFz == 0 or self.dim < 3: return + if self.nFz == 0 or self.dim < 3: + return return self._getTensorGrid('Fz') @property def gridEx(self): """Edge staggered grid in the x direction.""" - if self.nEx == 0: return + if self.nEx == 0: + return return self._getTensorGrid('Ex') @property def gridEy(self): """Edge staggered grid in the y direction.""" - if self.nEy == 0 or self.dim < 2: return + if self.nEy == 0 or self.dim < 2: + return return self._getTensorGrid('Ey') @property def gridEz(self): """Edge staggered grid in the z direction.""" - if self.nEz == 0 or self.dim < 3: return + if self.nEz == 0 or self.dim < 3: + return return self._getTensorGrid('Ez') def _getTensorGrid(self, key): @@ -220,22 +236,22 @@ def getTensor(self, key): """ - if key == 'Fx': - ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + if key == 'Fx': + ten = [self.vectorNx, self.vectorCCy, self.vectorCCz] elif key == 'Fy': - ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + ten = [self.vectorCCx, self.vectorNy, self.vectorCCz] elif key == 'Fz': - ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz] elif key == 'Ex': - ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + ten = [self.vectorCCx, self.vectorNy, self.vectorNz] elif key == 'Ey': - ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + ten = [self.vectorNx, self.vectorCCy, self.vectorNz] elif key == 'Ez': - ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + ten = [self.vectorNx, self.vectorNy, self.vectorCCz] elif key == 'CC': ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] elif key == 'N': - ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + ten = [self.vectorNx, self.vectorNy, self.vectorNz] return [t for t in ten if t is not None] @@ -254,14 +270,19 @@ def isInside(self, pts, locType='N'): tensors = self.getTensor(locType) if locType == 'N' and self._meshType == 'CYL': - #NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction! + # NOTE: for a CYL mesh we add a node to check if we are inside in + # the radial direction! tensors[0] = np.r_[0., tensors[0]] tensors[1] = np.r_[tensors[1], 2.0*np.pi] inside = np.ones(pts.shape[0], dtype=bool) for i, tensor in enumerate(tensors): TOL = np.diff(tensor).min() * 1.0e-10 - inside = inside & (pts[:, i] >= tensor.min()-TOL) & (pts[:, i] <= tensor.max()+TOL) + inside = ( + inside & + (pts[:, i] >= tensor.min()-TOL) & + (pts[:, i] <= tensor.max()+TOL) + ) return inside def _getInterpolationMat(self, loc, locType='CC', zerosOutside=False): @@ -293,7 +314,9 @@ def _getInterpolationMat(self, loc, locType='CC', zerosOutside=False): assert np.all(self.isInside(loc)), "Points outside of mesh" else: indZeros = np.logical_not(self.isInside(loc)) - loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')]) + loc[indZeros, :] = np.array([ + v.mean() for v in self.getTensor('CC') + ]) if locType in ['Fx', 'Fy', 'Fz', 'Ex', 'Ey', 'Ez']: ind = {'x': 0, 'y': 1, 'z': 2}[locType[1]] @@ -320,7 +343,8 @@ def _getInterpolationMat(self, loc, locType='CC', zerosOutside=False): else: raise NotImplementedError( - 'getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim) + 'getInterpolationMat: locType==' + locType + + ' and mesh.dim==' + str(self.dim) ) if zerosOutside: @@ -352,7 +376,9 @@ def getInterpolationMat(self, loc, locType='CC', zerosOutside=False): """ return self._getInterpolationMat(loc, locType, zerosOutside) - def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): + def _fastInnerProduct( + self, projType, prop=None, invProp=False, invMat=False + ): """ Fast version of getFaceInnerProduct. This does not handle the case of a full tensor prop. @@ -365,8 +391,9 @@ def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): :rtype: scipy.sparse.csr_matrix :return: M, the inner product matrix (nF, nF) """ - assert projType in ['F', 'E'], ("projType must be 'F' for faces or 'E'" - " for edges") + assert projType in ['F', 'E'], ( + "projType must be 'F' for faces or 'E' for edges" + ) if prop is None: prop = np.ones(self.nC) @@ -398,7 +425,7 @@ def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): # for faces, x, z matters, for edges, y (which is theta) matters if self._meshType == 'CYL': if projType == 'E': - prop = prop[:, 1] # this is the action of a projection mat + prop = prop[:, 1] # this is the action of a projection mat elif projType == 'F': prop = prop[:, [0, 2]] @@ -441,23 +468,26 @@ def _fastInnerProductDeriv(self, projType, prop, invProp=False, else: n_elements = self.dim - if tensorType == 0: # isotropic, constant Av = getattr(self, 'ave'+projType+'2CC') V = utils.sdiag(self.vol) - ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), - np.zeros(self.nC))), - shape=(self.nC, 1)) + ones = sp.csr_matrix( + (np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), + shape=(self.nC, 1) + ) if not invMat and not invProp: dMdprop = n_elements * Av.T * V * ones elif invMat and invProp: - dMdprop = n_elements * (utils.sdiag(MI.diagonal()**2) * Av.T * - V * ones * utils.sdiag(1./prop**2)) + dMdprop = n_elements * ( + utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * + utils.sdiag(1./prop**2) + ) elif invProp: dMdprop = n_elements * Av.T * V * utils.sdiag(- 1./prop**2) elif invMat: - dMdprop = n_elements * (utils.sdiag(- MI.diagonal()**2) * Av.T - * V) + dMdprop = n_elements * ( + utils.sdiag(- MI.diagonal()**2) * Av.T * V + ) elif tensorType == 1: # isotropic, variable in space Av = getattr(self, 'ave'+projType+'2CC') @@ -465,13 +495,16 @@ def _fastInnerProductDeriv(self, projType, prop, invProp=False, if not invMat and not invProp: dMdprop = n_elements * Av.T * V elif invMat and invProp: - dMdprop = n_elements * (utils.sdiag(MI.diagonal()**2) * Av.T * - V * utils.sdiag(1./prop**2)) + dMdprop = n_elements * ( + utils.sdiag(MI.diagonal()**2) * Av.T * V * + utils.sdiag(1./prop**2) + ) elif invProp: dMdprop = n_elements * Av.T * V * utils.sdiag(-1./prop**2) elif invMat: - dMdprop = n_elements * (utils.sdiag(- MI.diagonal()**2) * Av.T - * V) + dMdprop = n_elements * ( + utils.sdiag(- MI.diagonal()**2) * Av.T * V + ) elif tensorType == 2: # anisotropic Av = getattr(self, 'ave'+projType+'2CCV') @@ -514,8 +547,10 @@ def innerProductDeriv(v=None): return None -class TensorMesh(BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, - InnerProducts, TensorMeshIO): +class TensorMesh( + BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, + InnerProducts, TensorMeshIO +): """ TensorMesh is a mesh class that deals with tensor product meshes. @@ -605,7 +640,6 @@ def printH(hx, outStr=''): return outStr - # --------------- Geometries --------------------- @property def vol(self): @@ -620,7 +654,9 @@ def vol(self): self._vol = utils.mkvc(np.outer(vh[0], vh[1])) elif self.dim == 3: # Cell sizes in each direction - self._vol = utils.mkvc(np.outer(utils.mkvc(np.outer(vh[0], vh[1])), vh[2])) + self._vol = utils.mkvc( + np.outer(utils.mkvc(np.outer(vh[0], vh[1])), vh[2]) + ) return self._vol @property @@ -639,10 +675,18 @@ def area(self): area2 = np.outer(vh[0], np.ones(n[1]+1)) self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)] elif(self.dim == 3): - area1 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], vh[2]))) - area2 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - area3 = np.outer(vh[0], utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)] + area1 = np.outer( + np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], vh[2])) + ) + area2 = np.outer( + vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])) + ) + area3 = np.outer( + vh[0], utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))) + ) + self._area = np.r_[ + utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3) + ] return self._area @property @@ -661,10 +705,21 @@ def edge(self): l2 = np.outer(np.ones(n[0]+1), vh[1]) self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2)] elif(self.dim == 3): - l1 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) - l2 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - l3 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2), utils.mkvc(l3)] + l1 = np.outer( + vh[0], + utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))) + ) + l2 = np.outer( + np.ones(n[0]+1), + utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))) + ) + l3 = np.outer( + np.ones(n[0]+1), + utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])) + ) + self._edge = np.r_[ + utils.mkvc(l1), utils.mkvc(l2), utils.mkvc(l3) + ] return self._edge @property @@ -672,23 +727,23 @@ def faceBoundaryInd(self): """ Find indices of boundary faces in each direction """ - if self.dim==1: - indxd = (self.gridFx==min(self.gridFx)) - indxu = (self.gridFx==max(self.gridFx)) + if self.dim == 1: + indxd = (self.gridFx == min(self.gridFx)) + indxu = (self.gridFx == max(self.gridFx)) return indxd, indxu - elif self.dim==2: - indxd = (self.gridFx[:, 0]==min(self.gridFx[:, 0])) - indxu = (self.gridFx[:, 0]==max(self.gridFx[:, 0])) - indyd = (self.gridFy[:, 1]==min(self.gridFy[:, 1])) - indyu = (self.gridFy[:, 1]==max(self.gridFy[:, 1])) + elif self.dim == 2: + indxd = (self.gridFx[:, 0] == min(self.gridFx[:, 0])) + indxu = (self.gridFx[:, 0] == max(self.gridFx[:, 0])) + indyd = (self.gridFy[:, 1] == min(self.gridFy[:, 1])) + indyu = (self.gridFy[:, 1] == max(self.gridFy[:, 1])) return indxd, indxu, indyd, indyu - elif self.dim==3: - indxd = (self.gridFx[:, 0]==min(self.gridFx[:, 0])) - indxu = (self.gridFx[:, 0]==max(self.gridFx[:, 0])) - indyd = (self.gridFy[:, 1]==min(self.gridFy[:, 1])) - indyu = (self.gridFy[:, 1]==max(self.gridFy[:, 1])) - indzd = (self.gridFz[:, 2]==min(self.gridFz[:, 2])) - indzu = (self.gridFz[:, 2]==max(self.gridFz[:, 2])) + elif self.dim == 3: + indxd = (self.gridFx[:, 0] == min(self.gridFx[:, 0])) + indxu = (self.gridFx[:, 0] == max(self.gridFx[:, 0])) + indyd = (self.gridFy[:, 1] == min(self.gridFy[:, 1])) + indyu = (self.gridFy[:, 1] == max(self.gridFy[:, 1])) + indzd = (self.gridFz[:, 2] == min(self.gridFz[:, 2])) + indzu = (self.gridFz[:, 2] == max(self.gridFz[:, 2])) return indxd, indxu, indyd, indyu, indzd, indzu @property @@ -696,21 +751,21 @@ def cellBoundaryInd(self): """ Find indices of boundary faces in each direction """ - if self.dim==1: - indxd = (self.gridCC==min(self.gridCC)) - indxu = (self.gridCC==max(self.gridCC)) + if self.dim == 1: + indxd = (self.gridCC == min(self.gridCC)) + indxu = (self.gridCC == max(self.gridCC)) return indxd, indxu - elif self.dim==2: - indxd = (self.gridCC[:, 0]==min(self.gridCC[:, 0])) - indxu = (self.gridCC[:, 0]==max(self.gridCC[:, 0])) - indyd = (self.gridCC[:, 1]==min(self.gridCC[:, 1])) - indyu = (self.gridCC[:, 1]==max(self.gridCC[:, 1])) + elif self.dim == 2: + indxd = (self.gridCC[:, 0] == min(self.gridCC[:, 0])) + indxu = (self.gridCC[:, 0] == max(self.gridCC[:, 0])) + indyd = (self.gridCC[:, 1] == min(self.gridCC[:, 1])) + indyu = (self.gridCC[:, 1] == max(self.gridCC[:, 1])) return indxd, indxu, indyd, indyu - elif self.dim==3: - indxd = (self.gridCC[:, 0]==min(self.gridCC[:, 0])) - indxu = (self.gridCC[:, 0]==max(self.gridCC[:, 0])) - indyd = (self.gridCC[:, 1]==min(self.gridCC[:, 1])) - indyu = (self.gridCC[:, 1]==max(self.gridCC[:, 1])) - indzd = (self.gridCC[:, 2]==min(self.gridCC[:, 2])) - indzu = (self.gridCC[:, 2]==max(self.gridCC[:, 2])) + elif self.dim == 3: + indxd = (self.gridCC[:, 0] == min(self.gridCC[:, 0])) + indxu = (self.gridCC[:, 0] == max(self.gridCC[:, 0])) + indyd = (self.gridCC[:, 1] == min(self.gridCC[:, 1])) + indyu = (self.gridCC[:, 1] == max(self.gridCC[:, 1])) + indzd = (self.gridCC[:, 2] == min(self.gridCC[:, 2])) + indzu = (self.gridCC[:, 2] == max(self.gridCC[:, 2])) return indxd, indxu, indyd, indyu, indzd, indzu From d76cd60a216b940e11708eedf2e9e1a8dcdbc7b6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 5 Sep 2017 13:00:58 -0700 Subject: [PATCH 25/34] remove duplication in args / kwargs of h --- discretize/TensorMesh.py | 52 ++++++++++++++++------------------------ 1 file changed, 21 insertions(+), 31 deletions(-) diff --git a/discretize/TensorMesh.py b/discretize/TensorMesh.py index 2cf856cd8..c48d10c9e 100644 --- a/discretize/TensorMesh.py +++ b/discretize/TensorMesh.py @@ -36,37 +36,27 @@ def __init__(self, h=None, x0=None, **kwargs): h_in = h x0_in = x0 - # Cell widths - if h_in is not None: - # Sanity Checks - assert 'h' not in kwargs.keys(), ( - 'Only 1 descriptor of cell sizes can be used to instantiate ' - 'the class' + # Sanity Checks + assert type(h_in) in [list, tuple], 'h_in must be a list' + assert len(h_in) in [1, 2, 3], ( + 'h_in must be of dimension 1, 2, or 3' + ) + + # build h + h = list(range(len(h_in))) + for i, h_i in enumerate(h_in): + if utils.isScalar(h_i) and type(h_i) is not np.ndarray: + # This gives you something over the unit cube. + h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) + elif type(h_i) is list: + h_i = utils.meshTensor(h_i) + assert isinstance(h_i, np.ndarray), ( + "h[{0:d}] is not a numpy array.".format(i) ) - assert type(h_in) in [list, tuple], 'h_in must be a list' - assert len(h_in) in [1, 2, 3], ( - 'h_in must be of dimension 1, 2, or 3' + assert len(h_i.shape) == 1, ( + "h[{0:d}] must be a 1D numpy array.".format(i) ) - - # build h - h = list(range(len(h_in))) - for i, h_i in enumerate(h_in): - if utils.isScalar(h_i) and type(h_i) is not np.ndarray: - # This gives you something over the unit cube. - h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) - elif type(h_i) is list: - h_i = utils.meshTensor(h_i) - assert isinstance(h_i, np.ndarray), ( - "h[{0:d}] is not a numpy array.".format(i) - ) - assert len(h_i.shape) == 1, ( - "h[{0:d}] must be a 1D numpy array.".format(i) - ) - h[i] = h_i[:] # make a copy. - - else: - assert 'h' in kwargs.keys(), 'cell widths must be provided' - h = kwargs.pop('h') + h[i] = h_i[:] # make a copy. # Origin of the mesh x0 = np.zeros(len(h)) @@ -94,8 +84,8 @@ def __init__(self, h=None, x0=None, **kwargs): assert (n == np.array([x.size for x in h])).all(), ( "Dimension mismatch. The provided n doesn't " ) - - n = np.array([x.size for x in h]) + else: + n = np.array([x.size for x in h]) super(BaseTensorMesh, self).__init__( n, x0=x0 From 018d5850ccd9939a817085e2c98e2da8e1d8f7c6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 5 Sep 2017 14:49:20 -0700 Subject: [PATCH 26/34] more robust validation on n in the BaseMesh, check some end-cases in the testing to make sure errors are raised if the user messes with n. --- discretize/BaseMesh.py | 45 ++++++++++++++++++++++++++++++++++- discretize/TensorMesh.py | 2 +- discretize/TreeMesh.py | 4 ++-- tests/base/test_properties.py | 26 ++++++++++++++++++++ 4 files changed, 73 insertions(+), 4 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index c797b931c..7f9abf92e 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -38,6 +38,11 @@ def __init__(self, n, x0=None): # 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( @@ -45,8 +50,46 @@ def check_n_shape(self, change): "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_vs_n(self, change): + 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" + if len(self.n) != len(change['value']): raise Exception( "Dimension mismatch. x0 has length {} != len(n) which is " diff --git a/discretize/TensorMesh.py b/discretize/TensorMesh.py index c48d10c9e..20bb250d1 100644 --- a/discretize/TensorMesh.py +++ b/discretize/TensorMesh.py @@ -26,7 +26,7 @@ class BaseTensorMesh(BaseMesh): properties.Array( "widths of the tensor mesh in a single dimension", dtype=float, - shape=("*",) + shape=("*",), ), max_length=3 ) diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index be6b66580..b7d4e7cb6 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -110,7 +110,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): "cell number", min=0 ), - default=set() + default=set ) def __init__(self, h, x0=None, **kwargs): @@ -128,7 +128,7 @@ def __init__(self, h, x0=None, **kwargs): # self._levels = levels self._levelBits = int(np.ceil(np.sqrt(self.levels)))+1 - self.__dirty__ = True #: The numbering is dirty! + self.__dirty__ = True #: The numbering is dirty! if 'cells' in kwargs.keys(): self.cells = kwargs.pop('cells') diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index 6b795039e..eb5c717b6 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -73,6 +73,18 @@ def test_copy(self): compare_meshes(mesh0, mesh1) print('ok\n') + def test_base_updates(self): + with self.assertRaises(Exception): + self.mesh.n = None + + # check that if h has been set, we can't mess up n + with self.assertRaises(Exception): + self.mesh.n = [6, 5, 9] + + # can't change dimensionality of a mesh + with self.assertRaises(Exception): + self.mesh.n = [4, 5] + class CylTest(unittest.TestCase): @@ -131,6 +143,14 @@ def test_copy(self): compare_meshes(mesh0, mesh1) print('ok\n') + def test_base_updates(self): + with self.assertRaises(Exception): + self.mesh.n = None + + # check that if h has been set, we can't mess up n + with self.assertRaises(Exception): + self.mesh.n = [6, 5, 9] + class CurviTest(unittest.TestCase): @@ -160,7 +180,13 @@ def test_copy(self): compare_meshes(mesh0, mesh1) print('ok\n') + def test_base_updates(self): + with self.assertRaises(Exception): + self.mesh.n = None + # check that if h has been set, we can't mess up n + with self.assertRaises(Exception): + self.mesh.n = [6, 5, 9] if __name__ == '__main__': From 5060219363fa8c3d613b206c823bca6c74453cb5 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 5 Sep 2017 14:58:07 -0700 Subject: [PATCH 27/34] move to MeshIO. import load_mesh into base discretize namespace. deserialize onto BaseMesh rather than properties.HasProperties --- discretize/BaseMesh.py | 10 ++++++---- discretize/MeshIO.py | 22 +++++++++++++++++++++- discretize/__init__.py | 1 + discretize/utils/__init__.py | 10 +++++----- discretize/utils/meshutils.py | 23 +++-------------------- tests/base/test_properties.py | 8 ++++---- 6 files changed, 40 insertions(+), 34 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index 7f9abf92e..abfae88d3 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -2,7 +2,7 @@ import properties import os import json -from discretize import utils +from .utils.matutils import mkvc class BaseMesh(properties.HasProperties): @@ -11,6 +11,8 @@ class BaseMesh(properties.HasProperties): BaseMesh should be inherited by meshes with a regular structure. """ + _REGISTRY = {} + # Properties n = properties.Array( "number of cells in each direction (dim, )", @@ -703,7 +705,7 @@ def r(self, x, xType='CC', outType='CC', format='V'): 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" @@ -720,7 +722,7 @@ def outKernal(xx, nn): 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.""" @@ -733,7 +735,7 @@ def switchKernal(xx): 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 + 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] diff --git a/discretize/MeshIO.py b/discretize/MeshIO.py index 25ce4717a..f0d4538bd 100644 --- a/discretize/MeshIO.py +++ b/discretize/MeshIO.py @@ -1,8 +1,28 @@ import os +import json +import properties import numpy as np -from discretize import utils import six +from . import utils +from .BaseMesh import BaseMesh + + +def load_mesh(filename): + """ + Open a json file and load the mesh into the target class + + As long as there are no namespace conflicts, the target __class__ + will be stored on the properties.HasProperties registry and may be + fetched from there. + + :param str filename: name of file to read in + """ + with open(filename, 'r') as outfile: + jsondict = json.load(outfile) + data = BaseMesh.deserialize(jsondict, trusted=True) + return data + class TensorMeshIO(object): diff --git a/discretize/__init__.py b/discretize/__init__.py index b1880256d..9e6509962 100644 --- a/discretize/__init__.py +++ b/discretize/__init__.py @@ -3,6 +3,7 @@ from discretize.CylMesh import CylMesh from discretize.CurvilinearMesh import CurvilinearMesh from discretize import Tests +from discretize.MeshIO import load_mesh try: from discretize.TreeMesh import TreeMesh diff --git a/discretize/utils/__init__.py b/discretize/utils/__init__.py index d1c296b30..f400218b0 100644 --- a/discretize/utils/__init__.py +++ b/discretize/utils/__init__.py @@ -2,14 +2,14 @@ from .matutils import ( mkvc, sdiag, sdInv, speye, kron3, spzeros, ddx, av, - av_extrap, ndgrid, ind2sub, sub2ind, getSubArray, - inv3X3BlockDiagonal, inv2X2BlockDiagonal, TensorType, - makePropertyTensor, invPropertyTensor, Zero, - Identity + av_extrap, ndgrid, ind2sub, sub2ind, getSubArray, + inv3X3BlockDiagonal, inv2X2BlockDiagonal, TensorType, + makePropertyTensor, invPropertyTensor, Zero, + Identity ) from .codeutils import (isScalar, asArray_N_x_Dim) from .meshutils import ( - load_mesh, exampleLrmGrid, meshTensor, closestPoints, ExtractCoreMesh, + exampleLrmGrid, meshTensor, closestPoints, ExtractCoreMesh, random_model ) from .curvutils import volTetra, faceInfo, indexCube diff --git a/discretize/utils/meshutils.py b/discretize/utils/meshutils.py index c3b6691d0..6fb96dde5 100644 --- a/discretize/utils/meshutils.py +++ b/discretize/utils/meshutils.py @@ -1,11 +1,10 @@ import numpy as np -import properties -import json +import scipy.ndimage as ndi +import scipy.sparse as sp + from .matutils import ndgrid from .codeutils import asArray_N_x_Dim from .codeutils import isScalar -import scipy.ndimage as ndi -import scipy.sparse as sp import sys if sys.version_info < (3,): @@ -14,22 +13,6 @@ num_types = [int, float] -def load_mesh(filename): - """ - Open a json file and load the mesh into the target class - - As long as there are no namespace conflicts, the target __class__ - will be stored on the properties.HasProperties registry and may be - fetched from there. - - :param str filename: name of file to read in - """ - with open(filename, 'r') as outfile: - jsondict = json.load(outfile) - data = properties.HasProperties.deserialize(jsondict, trusted=True) - return data - - def exampleLrmGrid(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" diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index eb5c717b6..af30c93eb 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -61,7 +61,7 @@ def test_save_load(self): print('\nTesting save / load of Tensor Mesh ...') mesh0 = self.mesh f = mesh0.save() - mesh1 = discretize.utils.load_mesh(f) + mesh1 = discretize.load_mesh(f) compare_meshes(mesh0, mesh1) os.remove(f) print('ok\n') @@ -97,7 +97,7 @@ def test_save_load(self): print('\nTesting save / load of Cyl Mesh ...') mesh0 = self.mesh f = mesh0.save() - mesh1 = discretize.utils.load_mesh(f) + mesh1 = discretize.load_mesh(f) compare_meshes(mesh0, mesh1) os.remove(f) print('ok\n') @@ -131,7 +131,7 @@ def test_save_load(self): print('\nTesting save / load of Tree Mesh ...') mesh0 = self.mesh f = mesh0.save() - mesh1 = discretize.utils.load_mesh(f) + mesh1 = discretize.load_mesh(f) compare_meshes(mesh0, mesh1) os.remove(f) print('ok\n') @@ -168,7 +168,7 @@ def test_save_load(self): print('\nTesting save / load of Curvi Mesh ...') mesh0 = self.mesh f = mesh0.save() - mesh1 = discretize.utils.load_mesh(f) + mesh1 = discretize.load_mesh(f) compare_meshes(mesh0, mesh1) os.remove(f) print('ok\n') From 5980f4371d930c2db2dfd115488d563a27969d12 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 5 Sep 2017 17:01:26 -0700 Subject: [PATCH 28/34] update pypi credentials --- credentials.tar.gz.enc | Bin 304 -> 272 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/credentials.tar.gz.enc b/credentials.tar.gz.enc index 729b03663c5b5a9cdf17c25fc84034ece2f74ddf..19e87d25b17dfb96834044c765fdf17ce0a4d330 100644 GIT binary patch literal 272 zcmV+r0q_3)5gM{XVjVI8Wa=%uo|0zMDJ(s(#$Vi7CB&hWg@`by6TSSSYF{wDC{+t9 zlqg@0Y=%&Y-@}pP{Uc;N`m_i-D%pw21WjHmN9n>_o>ATcK{DDf&L zp^}d_?X)MQ-3jfnCz}GcZ>^UNDDnSvb`33G%Ta-Hs*VyogzIFFeszw1npfS?~ Wj?KM3e$^k`+CYoQ)wQQR0nW*->x$0+ literal 304 zcmV-00nh%i^@HbXFBZrk@IzwOD@GGJXd}AD9AufW=Fru|bo_`OO?m{^7<7@y8y)`b z5A|u&DO!`an8J!rQ$yZD0tjVu1<>``EF-vJ2h#6P<}U7Ol#{q$+woCLH>#T!M73dy`bTe$_e63lA3$sX2KLev8NN`0R1_{G2 zx86aiio;U03v Date: Tue, 5 Sep 2017 17:01:53 -0700 Subject: [PATCH 29/34] =?UTF-8?q?Bump=20version:=200.1.9=20=E2=86=92=200.1?= =?UTF-8?q?.10b0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- discretize/__init__.py | 2 +- docs/conf.py | 4 ++-- setup.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 84e9b99ec..f8a725304 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.9 +current_version = 0.1.10b0 files = setup.py discretize/__init__.py docs/conf.py diff --git a/discretize/__init__.py b/discretize/__init__.py index 9e6509962..02a56ebd1 100644 --- a/discretize/__init__.py +++ b/discretize/__init__.py @@ -18,7 +18,7 @@ """ ) -__version__ = '0.1.9' +__version__ = '0.1.10b0' __author__ = 'SimPEG Team' __license__ = 'MIT' __copyright__ = '2013 - 2017, SimPEG Developers, http://simpeg.xyz' diff --git a/docs/conf.py b/docs/conf.py index 9459bc10e..fcd26b469 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -60,9 +60,9 @@ # built documents. # # The short X.Y version. -version = '0.1.9' +version = '0.1.10b0' # The full version, including alpha/beta/rc tags. -release = '0.1.9' +release = '0.1.10b0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/setup.py b/setup.py index b9df2f9d0..00a27ef37 100644 --- a/setup.py +++ b/setup.py @@ -56,7 +56,7 @@ def configuration(parent_package='', top_path=None): setup( name="discretize", - version="0.1.9", + version="0.1.10b0", install_requires=[ 'numpy>=1.7', 'scipy>=0.13', From 23b5f45a8cf24e0985c7cf4d845d249d95960c5f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 10 Sep 2017 12:59:28 +0200 Subject: [PATCH 30/34] make n a private property --- discretize/BaseMesh.py | 36 +++++++++++++++++------------------ discretize/CurvilinearMesh.py | 4 ++-- tests/base/test_properties.py | 14 +++++++------- 3 files changed, 27 insertions(+), 27 deletions(-) diff --git a/discretize/BaseMesh.py b/discretize/BaseMesh.py index abfae88d3..3b3f1c311 100644 --- a/discretize/BaseMesh.py +++ b/discretize/BaseMesh.py @@ -14,7 +14,7 @@ class BaseMesh(properties.HasProperties): _REGISTRY = {} # Properties - n = properties.Array( + _n = properties.Array( "number of cells in each direction (dim, )", dtype=int, required=True, @@ -30,15 +30,15 @@ class BaseMesh(properties.HasProperties): # Instantiate the class def __init__(self, n, x0=None): - self.n = n # number of dimensions + self._n = n # number of dimensions if x0 is None: - self.x0 = np.zeros(len(self.n)) + self.x0 = np.zeros(len(self._n)) else: self.x0 = x0 # Validators - @properties.validator('n') + @properties.validator('_n') def check_n_shape(self, change): assert ( not isinstance(change['value'], properties.utils.Sentinel) and @@ -92,7 +92,7 @@ def check_x0(self, change): change['value'] is not None ), "n must be set prior to setting x0" - if len(self.n) != len(change['value']): + if len(self._n) != len(change['value']): raise Exception( "Dimension mismatch. x0 has length {} != len(n) which is " "{}".format(len(x0), len(n)) @@ -105,7 +105,7 @@ def dim(self): :rtype: int :return: dim """ - return len(self.n) + return len(self._n) @property def nC(self): @@ -122,7 +122,7 @@ def nC(self): M = discretize.TensorMesh([np.ones(n) for n in [2,3]]) M.plotGrid(centers=True, showIt=True) """ - return int(self.n.prod()) + return int(self._n.prod()) @property def nN(self): @@ -139,7 +139,7 @@ def nN(self): M = discretize.TensorMesh([np.ones(n) for n in [2,3]]) M.plotGrid(nodes=True, showIt=True) """ - return int((self.n+1).prod()) + return int((self._n+1).prod()) @property def nEx(self): @@ -148,7 +148,7 @@ def nEx(self): :rtype: int :return: nEx """ - return int((self.n + np.r_[0, 1, 1][:self.dim]).prod()) + return int((self._n + np.r_[0, 1, 1][:self.dim]).prod()) @property def nEy(self): @@ -159,7 +159,7 @@ def nEy(self): """ if self.dim < 2: return None - return int((self.n + np.r_[1, 0, 1][:self.dim]).prod()) + return int((self._n + np.r_[1, 0, 1][:self.dim]).prod()) @property def nEz(self): @@ -170,7 +170,7 @@ def nEz(self): """ if self.dim < 3: return None - return int((self.n + np.r_[1, 1, 0][:self.dim]).prod()) + return int((self._n + np.r_[1, 1, 0][:self.dim]).prod()) @property def vnE(self): @@ -209,7 +209,7 @@ def nFx(self): :rtype: int :return: nFx """ - return int((self.n + np.r_[1, 0, 0][:self.dim]).prod()) + return int((self._n + np.r_[1, 0, 0][:self.dim]).prod()) @property def nFy(self): @@ -220,7 +220,7 @@ def nFy(self): """ if self.dim < 2: return None - return int((self.n + np.r_[0, 1, 0][:self.dim]).prod()) + return int((self._n + np.r_[0, 1, 0][:self.dim]).prod()) @property def nFz(self): @@ -231,7 +231,7 @@ def nFz(self): """ if self.dim < 3: return None - return int((self.n + np.r_[0, 0, 1][:self.dim]).prod()) + return int((self._n + np.r_[0, 0, 1][:self.dim]).prod()) @property def vnF(self): @@ -386,7 +386,7 @@ def nCx(self): :rtype: int :return: nCx """ - return int(self.n[0]) + return int(self._n[0]) @property def nCy(self): @@ -397,7 +397,7 @@ def nCy(self): """ if self.dim < 2: return None - return int(self.n[1]) + return int(self._n[1]) @property def nCz(self): @@ -408,7 +408,7 @@ def nCz(self): """ if self.dim < 3: return None - return int(self.n[2]) + return int(self._n[2]) @property def vnC(self): @@ -727,7 +727,7 @@ def outKernal(xx, nn): def switchKernal(xx): """Switches over the different options.""" if xType in ['CC', 'N']: - nn = (self.n) if xType == 'CC' else (self.n+1) + nn = (self._n) if xType == 'CC' else (self._n+1) assert xx.size == np.prod(nn), ( "Number of elements must not change." ) diff --git a/discretize/CurvilinearMesh.py b/discretize/CurvilinearMesh.py index be1d7d5e4..2c802d1a2 100644 --- a/discretize/CurvilinearMesh.py +++ b/discretize/CurvilinearMesh.py @@ -62,8 +62,8 @@ def __init__(self, nodes=None, **kwargs): self.nodes = nodes - if 'n' in kwargs.keys(): - n = kwargs.pop('n') + if '_n' in kwargs.keys(): + n = kwargs.pop('_n') assert (n == np.array(self.nodes[0].shape)-1).all(), ( "Unexpected n-values. {} was provided, {} was expected".format( n, np.array(self.nodes[0].shape)-1 diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index af30c93eb..d22038397 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -75,15 +75,15 @@ def test_copy(self): def test_base_updates(self): with self.assertRaises(Exception): - self.mesh.n = None + self.mesh._n = None # check that if h has been set, we can't mess up n with self.assertRaises(Exception): - self.mesh.n = [6, 5, 9] + self.mesh._n = [6, 5, 9] # can't change dimensionality of a mesh with self.assertRaises(Exception): - self.mesh.n = [4, 5] + self.mesh._n = [4, 5] class CylTest(unittest.TestCase): @@ -145,11 +145,11 @@ def test_copy(self): def test_base_updates(self): with self.assertRaises(Exception): - self.mesh.n = None + self.mesh._n = None # check that if h has been set, we can't mess up n with self.assertRaises(Exception): - self.mesh.n = [6, 5, 9] + self.mesh._n = [6, 5, 9] class CurviTest(unittest.TestCase): @@ -182,11 +182,11 @@ def test_copy(self): def test_base_updates(self): with self.assertRaises(Exception): - self.mesh.n = None + self.mesh._n = None # check that if h has been set, we can't mess up n with self.assertRaises(Exception): - self.mesh.n = [6, 5, 9] + self.mesh._n = [6, 5, 9] if __name__ == '__main__': From 76d78c55c5fce8f89e46774c75d1f6d08291adc5 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 10 Sep 2017 13:10:52 +0200 Subject: [PATCH 31/34] cells is a private variable --- discretize/MeshIO.py | 2 +- discretize/TreeMesh.py | 44 +++++++++++++++++------------------ tests/base/test_properties.py | 4 ++-- 3 files changed, 25 insertions(+), 25 deletions(-) diff --git a/discretize/MeshIO.py b/discretize/MeshIO.py index f0d4538bd..f4b0c3a2b 100644 --- a/discretize/MeshIO.py +++ b/discretize/MeshIO.py @@ -485,7 +485,7 @@ def readUBC(TreeMesh, meshFile): # Make the tree mesh mesh = TreeMesh([h1, h2, h3], x0=x0) - mesh.cells = set([mesh._index(p) for p in simpegPointers.tolist()]) + mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) # Figure out the reordering mesh._simpegReorderUBC = np.argsort( diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index b7d4e7cb6..3f1b39ba8 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -104,7 +104,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): min=0 ) - cells = properties.Set( + _cells = properties.Set( "cells in the tree mesh", properties.Integer( "cell number", @@ -130,10 +130,10 @@ def __init__(self, h, x0=None, **kwargs): self.__dirty__ = True #: The numbering is dirty! - if 'cells' in kwargs.keys(): - self.cells = kwargs.pop('cells') + if '_cells' in kwargs.keys(): + self._cells = kwargs.pop('_cells') else: - self.cells.add(0) + self._cells.add(0) @properties.validator('levels') def check_levels(self, change): @@ -184,7 +184,7 @@ def fill(self): def maxLevel(self): """The maximum level used, which may be less than `levels`.""" l = 0 - for cell in self.cells: + for cell in self._cells: p = self._pointer(cell) l = max(l, p[-1]) return l @@ -232,7 +232,7 @@ def printH(hx, outStr=''): @property def nC(self): - return len(self.cells) + return len(self._cells) @property def nN(self): @@ -398,7 +398,7 @@ def ntEz(self): @property def _sortedCells(self): if getattr(self, '__sortedCells', None) is None: - self.__sortedCells = sorted(self.cells) + self.__sortedCells = sorted(self._cells) return self.__sortedCells @property @@ -439,7 +439,7 @@ def _pointer(self, index): return TreeUtils.point(self.dim, MAX_BITS, self._levelBits, index) def __contains__(self, v): - return self._asIndex(v) in self.cells + return self._asIndex(v) in self._cells def refine( self, function=None, recursive=True, cells=None, balance=True, @@ -454,7 +454,7 @@ def refine( self.__dirty__ = True if verbose: print('Refining Mesh') - cells = cells if cells is not None else sorted(self.cells) + cells = cells if cells is not None else sorted(self._cells) recurse = [] tic = time.time() for cell in cells: @@ -498,11 +498,11 @@ def corsen( self.__dirty__ = True if verbose: print('Corsening Mesh') - cells = cells if cells is not None else sorted(self.cells) + cells = cells if cells is not None else sorted(self._cells) recurse = [] tic = time.time() for cell in cells: - if cell not in self.cells: continue # already removed + if cell not in self._cells: continue # already removed p = self._pointer(cell) if p[-1] >= self.levels: continue result = function(Cell(self, cell, p)) @@ -537,8 +537,8 @@ def _refineCell(self, ind, pointer=None): raise IndexError(ind) children = self._childPointers(pointer, returnAll=True) for child in children: - self.cells.add(self._asIndex(child)) - self.cells.remove(ind) + self._cells.add(self._asIndex(child)) + self._cells.remove(ind) return [self._asIndex(child) for child in children] def _corsenCell(self, ind, pointer=None): @@ -549,9 +549,9 @@ def _corsenCell(self, ind, pointer=None): parent = self._parentPointer(pointer) children = self._childPointers(parent, returnAll=True) for child in children: - self.cells.remove(self._asIndex(child)) + self._cells.remove(self._asIndex(child)) parentInd = self._asIndex(parent) - self.cells.add(parentInd) + self._cells.add(parentInd) return [parentInd] def _asPointer(self, ind): @@ -676,7 +676,7 @@ def balance(self, recursive=True, cells=None, verbose=False, _inRecursion=False) self.__dirty__ = True if verbose: print('Balancing Mesh:') - cells = cells if cells is not None else sorted(self.cells) + cells = cells if cells is not None else sorted(self._cells) # calcDepth = lambda i: lambda A: i if type(A) is not list else max(map(calcDepth(i+1), A)) # flatten = lambda A: A if calcDepth(0)(A) == 1 else flatten([_ for __ in A for _ in (__ if type(__) is list else [__])]) @@ -716,7 +716,7 @@ def balance(self, recursive=True, cells=None, verbose=False, _inRecursion=False) @property def gridCC(self): if getattr(self, '_gridCC', None) is None: - self._gridCC = np.zeros((len(self.cells),self.dim)) + self._gridCC = np.zeros((len(self._cells),self.dim)) for ii, ind in enumerate(self._sortedCells): p = self._asPointer(ind) self._gridCC[ii, :] = self._cellC(p) + self.x0 @@ -775,7 +775,7 @@ def gridEz(self): @property def vol(self): if getattr(self, '_vol', None) is None: - self._vol = np.zeros(len(self.cells)) + self._vol = np.zeros(len(self._cells)) for ii, ind in enumerate(self._sortedCells): p = self._asPointer(ind) self._vol[ii] = np.prod(self._cellH(p)) @@ -818,7 +818,7 @@ def _createNumberingSets(self, force=False): self._edgesZ = set() - for ind in self.cells: + for ind in self._cells: p = self._asPointer(ind) w = self._levelWidth(p[-1]) if self.dim == 2: @@ -889,7 +889,7 @@ def _numberCells(self, force=False): return self._cc2i = dict() self._i2cc = dict() - for ii, c in enumerate(sorted(self.cells)): + for ii, c in enumerate(sorted(self._cells)): self._cc2i[c] = ii self._i2cc[ii] = c self.__dirtyCells__ = False @@ -912,7 +912,7 @@ def _numberFaces(self, force=False): return self._createNumberingSets(force=force) - for ind in self.cells: + for ind in self._cells: p = self._asPointer(ind) w = self._levelWidth(p[-1]) @@ -2284,7 +2284,7 @@ def doSlice(v): _ = antiNormalInd X = [] Y = [] - for cell in self.cells: + for cell in self._cells: p = self._pointer(cell) n, h = self._cellN(p), self._cellH(p) if n[normalInd]indLoc: diff --git a/tests/base/test_properties.py b/tests/base/test_properties.py index d22038397..9f987811b 100644 --- a/tests/base/test_properties.py +++ b/tests/base/test_properties.py @@ -40,8 +40,8 @@ def compare_meshes(mesh0, mesh1): assert (mesh0.vol == mesh1.vol).all() # Tree mesh specific - if hasattr(mesh0, 'cells'): - assert (mesh0.cells == mesh1.cells) + if hasattr(mesh0, '_cells'): + assert (mesh0._cells == mesh1._cells) # curvi-specific if hasattr(mesh0, 'nodes'): From d15dda39789af92a43f6662e1dfa3d483e6d919b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 10 Sep 2017 13:12:33 +0200 Subject: [PATCH 32/34] update requirements to rely on properties 0.3.6b0 --- requirements_dev.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements_dev.txt b/requirements_dev.txt index 5eef8a0bc..aacca4af0 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -12,6 +12,6 @@ cython pymatsolver>=0.1.2 ipython matplotlib -properties[math] +properties[math]>=0.3.6b0 wheel twine From ad2d77adf0146b33993a47f977b0ad853b68f6c6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 10 Sep 2017 13:20:53 +0200 Subject: [PATCH 33/34] levels is a private property of the tree mesh --- discretize/TreeMesh.py | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/discretize/TreeMesh.py b/discretize/TreeMesh.py index 3f1b39ba8..4f1fb757e 100644 --- a/discretize/TreeMesh.py +++ b/discretize/TreeMesh.py @@ -99,7 +99,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): _meshType = 'TREE' - levels = properties.Integer( + _levels = properties.Integer( "discretization level", min=0 ) @@ -117,16 +117,16 @@ def __init__(self, h, x0=None, **kwargs): assert type(h) is list, 'h must be a list' assert len(h) in [2, 3], "TreeMesh is only in 2D or 3D." - if 'levels' in kwargs.keys(): - self.levels = kwargs.pop('levels') + if '_levels' in kwargs.keys(): + self._levels = kwargs.pop('_levels') BaseTensorMesh.__init__(self, h, x0, **kwargs) - if self.levels is None: - self.levels = int(np.log2(len(self.h[0]))) + if self._levels is None: + self._levels = int(np.log2(len(self.h[0]))) # self._levels = levels - self._levelBits = int(np.ceil(np.sqrt(self.levels)))+1 + self._levelBits = int(np.ceil(np.sqrt(self._levels)))+1 self.__dirty__ = True #: The numbering is dirty! @@ -135,7 +135,7 @@ def __init__(self, h, x0=None, **kwargs): else: self._cells.add(0) - @properties.validator('levels') + @properties.validator('_levels') def check_levels(self, change): assert np.all( len(_) == 2**change['value'] for _ in self.h @@ -431,7 +431,7 @@ def permuteE(self): def _index(self, pointer): assert len(pointer) is self.dim+1 - assert pointer[-1] <= self.levels + assert pointer[-1] <= self._levels return TreeUtils.index(self.dim, MAX_BITS, self._levelBits, pointer[:-1], pointer[-1]) def _pointer(self, index): @@ -459,7 +459,7 @@ def refine( tic = time.time() for cell in cells: p = self._pointer(cell) - if p[-1] >= self.levels: continue + if p[-1] >= self._levels: continue result = function(Cell(self, cell, p)) if type(result) is bool: do = result @@ -504,7 +504,7 @@ def corsen( for cell in cells: if cell not in self._cells: continue # already removed p = self._pointer(cell) - if p[-1] >= self.levels: continue + if p[-1] >= self._levels: continue result = function(Cell(self, cell, p)) if type(result) is bool: do = result @@ -559,7 +559,7 @@ def _asPointer(self, ind): return self._pointer(ind) if type(ind) is list: assert len(ind) == (self.dim + 1), str(ind) +' is not valid pointer' - assert ind[-1] <= self.levels, str(ind) +' is not valid pointer' + assert ind[-1] <= self._levels, str(ind) +' is not valid pointer' return ind if isinstance(ind, np.ndarray): return ind.tolist() @@ -626,12 +626,12 @@ def _cellC(self, p): return (np.array(self._cellH(p))/2.0 + self._cellN(p)).tolist() def _levelWidth(self, level): - return 2**(self.levels - level) + return 2**(self._levels - level) def _isInsideMesh(self, pointer): inside = True for p in pointer[:-1]: - inside = inside and p >= 0 and p < 2**self.levels + inside = inside and p >= 0 and p < 2**self._levels return inside def _getNextCell(self, ind, direction=0, positive=True, _lookUp=True): @@ -643,7 +643,7 @@ def _getNextCell(self, ind, direction=0, positive=True, _lookUp=True): if direction >= self.dim: return None pointer = self._asPointer(ind) - if pointer[-1] > self.levels: + if pointer[-1] > self._levels: return None step = (1 if positive else -1) * self._levelWidth(pointer[-1]) @@ -656,7 +656,7 @@ def _getNextCell(self, ind, direction=0, positive=True, _lookUp=True): if nextCell in self: return self._index(nextCell) - if nextCell[-1] + 1 <= self.levels: # if I am not the smallest. + if nextCell[-1] + 1 <= self._levels: # if I am not the smallest. children = self._childPointers(pointer, direction=direction, positive=positive) nextCells = [self._getNextCell(child, direction=direction, positive=positive, _lookUp=False) for child in children] if nextCells[0] is not None: @@ -685,7 +685,7 @@ def balance(self, recursive=True, cells=None, verbose=False, _inRecursion=False) for cell in cells: p = self._asPointer(cell) - if p[-1] == self.levels: continue + if p[-1] == self._levels: continue cs = list(range(6)) cs[0] = self._getNextCell(cell, direction=0, positive=False) @@ -1033,7 +1033,7 @@ def _hanging(self, force=False): # Compute from x faces for fx in self._facesX: p = self._pointer(fx) - if p[-1] + 1 > self.levels: continue + if p[-1] + 1 > self._levels: continue sl = p[-1] + 1 #: small level test = self._index(p[:-1] + [sl]) if test not in self._facesX: @@ -1130,7 +1130,7 @@ def _hanging(self, force=False): # Compute from y faces for fy in self._facesY: p = self._pointer(fy) - if p[-1] + 1 > self.levels: continue + if p[-1] + 1 > self._levels: continue sl = p[-1] + 1 #: small level test = self._index(p[:-1] + [sl]) if test not in self._facesY: @@ -1230,7 +1230,7 @@ def _hanging(self, force=False): # Compute from z faces for fz in self._facesZ: p = self._pointer(fz) - if p[-1] + 1 > self.levels: continue + if p[-1] + 1 > self._levels: continue sl = p[-1] + 1 #: small level test = self._index(p[:-1] + [sl]) if test not in self._facesZ: @@ -1924,7 +1924,7 @@ def point2index(self, locs): out = [] for pointer in zip(*pointers): - for level in range(self.levels+1): + for level in range(self._levels+1): width = self._levelWidth(level) testPointer = [((p-1)//width)*width for p in pointer] + [level] test = self._index(testPointer) From a75afa59d19092c7b03f0f9df86f1f9f5319ced4 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 11 Sep 2017 09:29:03 +0200 Subject: [PATCH 34/34] use the array property with shapes as a set for nodes in the curvi mesh. See updates in aranzgeo/properties#195 --- discretize/CurvilinearMesh.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/discretize/CurvilinearMesh.py b/discretize/CurvilinearMesh.py index 2c802d1a2..d829d35e8 100644 --- a/discretize/CurvilinearMesh.py +++ b/discretize/CurvilinearMesh.py @@ -47,12 +47,9 @@ class CurvilinearMesh( nodes = properties.List( "List of arrays describing the node locations", - prop=properties.Union( + prop=properties.Array( "node locations in an n-dimensional array", - props=[ - properties.Array('2D array', shape=('*', '*')), - properties.Array('3D array', shape=('*', '*', '*')), - ] + shape={('*', '*'), ('*', '*', '*')} ), min_length=2, max_length=3