From 13325e63fc3b986c228b03c62da9d52317e794ef Mon Sep 17 00:00:00 2001 From: Luz Angelica Caudillo Mata Date: Tue, 2 Jul 2013 12:59:02 -0700 Subject: [PATCH 1/6] Implements the first draft of the TensorMesh class with basic visualization functions. It works for 1, 2 and 3 dimensions. Provides nodal, cell centers and staggered grids. --- code/TensorMesh.py | 257 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 257 insertions(+) create mode 100644 code/TensorMesh.py diff --git a/code/TensorMesh.py b/code/TensorMesh.py new file mode 100644 index 00000000..3a061561 --- /dev/null +++ b/code/TensorMesh.py @@ -0,0 +1,257 @@ +#------------------------------------------------------------------------------- +# Packages +#------------------------------------------------------------------------------- +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + +#------------------------------------------------------------------------------- +# Class definition +#------------------------------------------------------------------------------- +class TensorMesh: + """ + Define nodal, cell-centered and staggered tensor meshes for 1, 2 and 3 + dimensions. + """ + + # ---------------------- Properties -------------------------------------- + h = None # Array Spacing or cell-sizes in each direction + x0 = None # Array Origin (x1,x2,x3) + dim = None # Int Dimension + n = None # Array Number of cells in each direction + nC = None # Int Total number of cells + nE = None # Int Total number of edges + nF = None # Int Total number of faces + + # ----------------------- Methods ---------------------------------------- + def __init__(self,h,x0): + """Compute number of edges,faces and cell-centers """ + + # Assign values to properties + self.h = h + self.x0 = x0 + + + #Compute derived properties + self.dim = np.size(x0) + + #Compute the num of cells in each direction + self.n = np.zeros((self.dim,1)) + + for d in range(self.dim): + self.n[d] = np.size(h[d]) + + # Compute the number of cell-centers + self.nC = np.prod(self.n) + + # Compute the number of edges (makes sense only for 3D) + # Equivalent to: + # nEdges = n[0] * (n[1]+1) * (n[2]+1) + # + (n[0]+1) * ny[1] * (nz[2]+1) + # + (n[0]+1) * (ny[1]+1)* (nz[2]) + + if self.dim == 3: + self.nE = np.prod(np.kron(np.ones((3,1)),self.n.T)+np.ones((3,3))-np.eye(3),1) + print self.nE + + # Compute the number of faces (makes sense only for 2 and 3D) + # Equivalent to + # nFaces = (n[0]+1) * n[1] * n[2] + # + n[0] * (ny[1]+1) * nz[2] + # + n[0] * ny[1] * (nz[2]+1) + if self.dim >=2: + self.nF = np.prod(np.kron(np.ones((self.dim,1)),self.n.T)+np.eye(self.dim),1) + print self.nF + + def xin(self,i): + """Construct the 1D nodal mesh from the ith-component of h. Return an array.""" + return np.insert(np.cumsum(self.h[i-1]),0,0.0) + self.x0[i-1] + + def xic(self,i): + """Construct the 1D cell-centerd mesh from the ith-component of h. Return an array.""" + return .5*( np.insert(np.cumsum(self.h[i-1][:,0:-1]),0,0.0) + np.cumsum(self.h[i-1])) + + def getNodalGrid(self): + """Construct nodal grid for 1, 2 and 3 dimensions""" + if self.dim==1: + return [self.xin(1)] + elif self.dim==2: + return self.ndgrid([self.xin(1),self.xin(2)]) + elif self.dim==3: + return self.ndgrid([self.xin(1),self.xin(2),self.xin(3)]) + + def ndgrid(self, xin): + """Form tensorial grid for 1, 2 and 3 dimensions. Return X1,X2,X3 arrays depending on the dimension""" + ei = lambda i : np.ones((np.size(xin[i-1]),1)) + + if self.dim==1: + return [xin] + elif self.dim==2: + X1 = np.kron(ei(2),xin[0]).reshape(-1,1) + X2 = np.kron(xin[1],ei(1).T).reshape(-1,1) + return X1,X2 + elif self.dim==3: + X1 = np.kron(ei(3),np.kron(ei(2),xin[0])).reshape(-1,1) + X2 = np.kron(ei(3).T,np.kron(xin[1],ei(1).T)).reshape(-1,1) + X3 = np.kron(xin[2],np.kron(ei(2),ei(1))).T.reshape(-1,1) + + return X1,X2,X3 + + def getCellCenteredGrid(self): + """Construct cell-centered grid for 1, 2 and 3 dimensions.""" + + if self.dim==1: + return [self.xic(1)] + elif self.dim==2: + return self.ndgrid([self.xic(1),self.xic(2)]) + elif self.dim==3: + return self.ndgrid([self.xic(1),self.xic(2),self.xic(3)]) + + def getFaceStgGrid(self,direction): + """Construct the face staggered grids for 2 and 3 dimensions.""" + + if self.dim==1: + print 'Error: dimension must be larger than 1' + elif self.dim==2: + if direction == 1: + return self.ndgrid([self.xin(1),self.xic(2)]) + elif direction == 2: + return self.ndgrid([self.xic(1),self.xin(2)]) + else: + print 'Error: direction must be equal to 1 or 2' + elif self.dim==3: + if direction == 1: + return self.ndgrid([self.xin(1),self.xic(2),self.xic(3)]) + elif direction == 2: + return self.ndgrid([self.xic(1),self.xin(2),self.xic(3)]) + elif direction == 3: + return self.ndgrid([self.xic(1),self.xic(2),self.xin(3)]) + else: + print 'Error: direction must be equal to 1, 2 or 3' + + def getEdgeStgGrid(self,direction): + """Construct the edge staggered grids for 3 dimension case.""" + if self.dim != 3: + print 'Error: dimension must be equal to 3' + else: + if direction == 1: + return self.ndgrid([self.xic(1),self.xin(2),self.xin(3)]) + elif direction == 2: + return self.ndgrid([self.xin(1),self.xic(2),self.xin(3)]) + elif direction == 3: + return self.ndgrid([self.xin(1),self.xin(2),self.xic(3)]) + else: + print 'Error: direction must be equal to 1, 2 or 3' + + def plotImage(self,I): + + if self.dim==1: + fig = plt.figure(1) + fig.clf() + ax=plt.subplot(111) + if np.size(I)==self.n[0]: + print 'cell-centered image' + xx = self.getCellCenteredGrid() + ax.plot(xx[0],I,'ro') + elif np.size(I)==self.n[0]+1: + print 'nodal image' + xx = self.getNodalGrid() + ax.plot(xx[0],I,'bs') + + fig.show() + + + def plotGrid(self): + """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.""" + if self.dim == 1: + fig = plt.figure(1) + fig.clf() + ax = plt.subplot(111) + xn = self.getNodalGrid() + xc = self.getCellCenteredGrid() + print xn + ax.hold(True) + ax.plot(xn,np.ones(np.shape(xn)),'bs') + ax.plot(xc,np.ones(np.shape(xc)),'ro') + ax.plot(xn,np.ones(np.shape(xn)),'k--') + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + fig.show() + elif self.dim == 2: + fig = plt.figure(2) + fig.clf() + ax = plt.subplot(111) + xn = self.getNodalGrid() + xc = self.getCellCenteredGrid() + xs1 = self.getFaceStgGrid(1) + xs2 = self.getFaceStgGrid(2) + + ax.hold(True) + ax.plot(xn[0],xn[1],'bs') + ax.plot(xc[0],xc[1],'ro') + ax.plot(xs1[0],xs1[1],'g>') + ax.plot(xs2[0],xs2[1],'g^') + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + fig.show() + elif self.dim == 3: + fig = plt.figure(3) + fig.clf() + ax = fig.add_subplot(111, projection='3d') + xn = self.getNodalGrid() + xc = self.getCellCenteredGrid() + xfs1 = self.getFaceStgGrid(1) + xfs2 = self.getFaceStgGrid(2) + xfs3 = self.getFaceStgGrid(3) + xes1 = self.getEdgeStgGrid(1) + xes2 = self.getEdgeStgGrid(2) + xes3 = self.getEdgeStgGrid(3) + + ax.hold(True) + ax.plot(xn[0],xn[1],'bs',zs=xn[2]) + ax.plot(xc[0],xc[1],'ro',zs=xc[2]) + ax.plot(xfs1[0],xfs1[1],'g>',zs=xfs1[2]) + ax.plot(xfs2[0],xfs2[1],'g<',zs=xfs2[2]) + ax.plot(xfs3[0],xfs3[1],'g^',zs=xfs3[2]) + ax.plot(xes1[0],xes1[1],'k>',zs=xes1[2]) + ax.plot(xes2[0],xes2[1],'k<',zs=xes2[2]) + ax.plot(xes3[0],xes3[1],'k^',zs=xes3[2]) + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + ax.set_zlabel('x3') + fig.show() + + + +if __name__ == '__main__': + print('Welcome to tensor mesh!') + + testDim = 1 + h1 = 0.3*np.ones((1,7)) + h1[:,0] = 0.5 + h1[:,-1] = 0.6 + h2 = .5* np.ones((1,4)) + h3 = .4* np.ones((1,6)) + x0 = np.zeros((3,1)) + + if testDim == 1: + h = [h1] + x0 = x0[0] + elif testDim==2: + h = [h1,h2] + x0 = x0[0:2] + else: + h = [h1,h2,h3] + + I = np.linspace(0,1,8) + M = TensorMesh(h,x0) + + xn = M.plotGrid() + + + \ No newline at end of file From e9d41ab95fd210562d6f23b9da02a9346172b1c5 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 5 Jul 2013 16:00:19 -0700 Subject: [PATCH 2/6] Changed Mesh to BaseMesh. Added conventions of n{C-F-E}{x-y-z} so mesh.nCx will return the number of cells in the x direction. Added a unitTest to ensure the mesh dimensions. --- code/BaseMesh.py | 179 ++++++++++++++++++++++++++++++++++++ code/Mesh.py | 130 -------------------------- code/tests/test_basemesh.py | 88 ++++++++++++++++++ code/tests/test_mesh.py | 74 --------------- 4 files changed, 267 insertions(+), 204 deletions(-) create mode 100644 code/BaseMesh.py delete mode 100644 code/Mesh.py create mode 100644 code/tests/test_basemesh.py delete mode 100644 code/tests/test_mesh.py diff --git a/code/BaseMesh.py b/code/BaseMesh.py new file mode 100644 index 00000000..90ad07f3 --- /dev/null +++ b/code/BaseMesh.py @@ -0,0 +1,179 @@ +import numpy as np + + +class BaseMesh(object): + """BaseMesh does all the counting you don't want to do.""" + def __init__(self, n, x0=None): + + # Check inputs + if x0 is None: + x0 = np.zeros(len(n)) + + if not len(n) == len(x0): + 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).ravel() + self._dim = len(n) + + def x0(): + doc = "Origin of the mesh" + fget = lambda self: self._x0 + return locals() + x0 = property(**x0()) + + def n(): + doc = "Number of Cells in each dimension (array of integers)" + fget = lambda self: self._n + return locals() + n = property(**n()) + + def dim(): + doc = "The dimension of the mesh (1, 2, or 3)." + fget = lambda self: self._dim + return locals() + dim = property(**dim()) + + def nCx(): + doc = "Number oc cells in the x direction" + fget = lambda self: self.n[0] + return locals() + nCx = property(**nCx()) + + def nCy(): + doc = "Number of cells in the y direction" + + def fget(self): + if self.dim > 1: + return self.n[1] + else: + return None + return locals() + nCy = property(**nCy()) + + def nCz(): + doc = "Number of cells in the z direction" + + def fget(self): + if self.dim > 2: + return self.n[2] + else: + return None + return locals() + nCz = property(**nCz()) + + def nC(): + doc = "Total number of cells" + fget = lambda self: np.prod(self.n) + return locals() + nC = property(**nC()) + + def nNx(): + doc = "Number of nodes in the x-direction" + fget = lambda self: self.nCx + 1 + return locals() + nNx = property(**nNx()) + + def nNy(): + doc = "Number of noes in the y-direction" + + def fget(self): + if self.dim > 1: + return self.n[1] + 1 + else: + return None + return locals() + nNy = property(**nNy()) + + def nNz(): + doc = "Number of nodes in the z-direction" + + def fget(self): + if self.dim > 2: + return self.n[2] + 1 + else: + return None + return locals() + nNz = property(**nNz()) + + def nN(): + doc = "Total number of nodes" + fget = lambda self: self.n + 1 + return locals() + nN = property(**nN()) + + def nEx(): + doc = "Number of x-edges" + fget = lambda self: np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None]) + return locals() + nEx = property(**nEx()) + + def nEy(): + doc = "Number of y-edges" + + def fget(self): + if self.dim > 1: + return np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None]) + else: + return None + return locals() + nEy = property(**nEy()) + + def nEz(): + doc = "Number of z-edges" + + def fget(self): + if self.dim > 2: + return np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None]) + else: + return None + return locals() + nEz = property(**nEz()) + + def nE(): + doc = "Total number of edges" + fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None]) + return locals() + nE = property(**nE()) + + def nFx(): + doc = "Number of x-faces" + fget = lambda self: np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None]) + return locals() + nFx = property(**nFx()) + + def nFy(): + doc = "Number of y-faces" + + def fget(self): + if self.dim > 1: + return np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None]) + else: + return None + return locals() + nFy = property(**nFy()) + + def nFz(): + doc = "Number of z-faces" + + def fget(self): + if self.dim > 2: + return np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None]) + else: + return None + return locals() + nFz = property(**nFz()) + + def nF(): + doc = "Total number of faces in each dimension" + fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None]) + return locals() + nF = property(**nF()) + +if __name__ == '__main__': + m = BaseMesh([3, 2, 4]) + print m.n diff --git a/code/Mesh.py b/code/Mesh.py deleted file mode 100644 index 8ca87165..00000000 --- a/code/Mesh.py +++ /dev/null @@ -1,130 +0,0 @@ -import numpy as np - - -class Mesh(object): - """docstring for Mesh""" - def __init__(self, h): - - if type(h) != list: - raise Exception("Type of h must be a list variable. e.g. [5, 4, 2] or [[1,1,1],[0.5,0.5]]") - - if np.sum([np.size(x) for x in h]) == len(h): - # We have specified a shorthand for the mesh e.g. [5, 4, 2] - # We will recreate the h, such that it lies on the unit cube/square/line - domain = 1. # (must be a float) - h = [np.ones(x)*(domain/x) for x in h] - - dim = len(h) - - if dim > 1 and np.all([len(np.shape(x)) > 1 and np.shape(x)[1] > 1 for x in h]): - # The h has internal structure, and is not a vector - # Thus, we must be describing the verticies of the mesh - # Hence, the mesh is a Logically Orthogonal Mesh - self.meshType = 'LOM' - else: - # Could add other checks, but here the default is a rectangular mesh - self.meshType = 'RECT' - - if self.meshType != 'LOM': - # Ensure that the h is a numpy array, with shape: (n,) - h = [np.array(x).ravel() for x in h] - - # Define the number of nodes - if self.meshType == 'LOM': - self._nnodes = np.array(np.shape(h[0])) - else: - self._nnodes = np.array([len(x) for x in h]) + 1 - - self._nc = self._nnodes - 1 - self._ncells = np.prod(self._nc) - self._h = h - self._dim = dim - - m = self._nnodes - if dim == 1: - self._nfaces = np.prod(m) - self._nedges = np.prod(m) - elif dim == 2: - self._nfx = m - [0, 1] - self._nfy = m - [1, 0] - self._nex = m - [1, 0] - self._ney = m - [0, 1] - - self._nfaces = [np.prod(self.nfx), np.prod(self.nfy)] - self._nedges = [np.prod(self.nex), np.prod(self.ney)] - elif dim == 3: - self._nfx = m - [0, 1, 1] - self._nfy = m - [1, 0, 1] - self._nfz = m - [1, 1, 0] - self._nex = m - [1, 0, 0] - self._ney = m - [0, 1, 0] - self._nez = m - [0, 0, 1] - - self._nfaces = [np.prod(self.nfx), np.prod(self.nfy), np.prod(self.nfz)] - self._nedges = [np.prod(self.nex), np.prod(self.ney), np.prod(self.nez)] - - def dim(): - doc = "The dimension of the mesh: 1, 2, or 3" - fget = lambda self: self._dim - return locals() - dim = property(**dim()) - - def nc(): - doc = "Number of cells in each direction of the mesh" - fget = lambda self: self._nc - return locals() - nc = property(**nc()) - - def ncells(): - doc = "Number of cells in the mesh" - fget = lambda self: self._ncells - return locals() - ncells = property(**ncells()) - - def nfaces(): - doc = "Number of faces in each direction of the mesh" - fget = lambda self: self._nfaces - return locals() - nfaces = property(**nfaces()) - - def nedges(): - doc = "Number of edges in each direction of the mesh" - fget = lambda self: self._nedges - return locals() - nedges = property(**nedges()) - - def nfx(): - doc = "Number of faces in the x direction of the mesh" - fget = lambda self: self._nfx if self.dim > 1 else None - return locals() - nfx = property(**nfx()) - - def nfy(): - doc = "Number of faces in the y direction of the mesh" - fget = lambda self: self._nfy if self.dim > 1 else None - return locals() - nfy = property(**nfy()) - - def nfz(): - doc = "Number of faces in the z direction of the mesh" - fget = lambda self: self._nfz if self.dim > 2 else None - return locals() - nfz = property(**nfz()) - - def nex(): - doc = "Number of edges in the x direction of the mesh" - fget = lambda self: self._nex if self.dim > 1 else None - return locals() - nex = property(**nex()) - - def ney(): - doc = "Number of edges in the y direction of the mesh" - fget = lambda self: self._ney if self.dim > 1 else None - return locals() - ney = property(**ney()) - - def nez(): - doc = "Number of edges in the z direction of the mesh" - fget = lambda self: self._nez if self.dim > 2 else None - return locals() - nez = property(**nez()) diff --git a/code/tests/test_basemesh.py b/code/tests/test_basemesh.py new file mode 100644 index 00000000..88655c60 --- /dev/null +++ b/code/tests/test_basemesh.py @@ -0,0 +1,88 @@ +import unittest +import sys +sys.path.append('../') +from BaseMesh import BaseMesh +import numpy as np + + +class TestMeshNumbers3D(unittest.TestCase): + + def setUp(self): + self.mesh = BaseMesh([6, 2, 3]) + + def test_meshDimensions(self): + self.assertTrue(self.mesh.dim, 3) + + def test_mesh_nc(self): + self.assertTrue(np.all(self.mesh.n == [6, 2, 3])) + + def test_mesh_nc_xyz(self): + x = np.all(self.mesh.nCx == 6) + y = np.all(self.mesh.nCy == 2) + z = np.all(self.mesh.nCz == 3) + + self.assertTrue(np.all([x, y, z])) + + def test_mesh_nf(self): + x = np.all(self.mesh.nFx == [7, 2, 3]) + y = np.all(self.mesh.nFy == [6, 3, 3]) + z = np.all(self.mesh.nFz == [6, 2, 4]) + + self.assertTrue(np.all([x, y, z])) + + def test_mesh_ne(self): + x = np.all(self.mesh.nEx == [6, 3, 4]) + y = np.all(self.mesh.nEy == [7, 2, 4]) + z = np.all(self.mesh.nEz == [7, 3, 3]) + + self.assertTrue(np.all([x, y, z])) + + def test_mesh_numbers(self): + c = self.mesh.nC == 36 + f = np.all(self.mesh.nF == [42, 54, 48]) + e = np.all(self.mesh.nE == [72, 56, 63]) + + self.assertTrue(np.all([c, f, e])) + + +class TestMeshNumbers2D(unittest.TestCase): + + def setUp(self): + self.mesh = BaseMesh([6, 2]) + + def test_meshDimensions(self): + self.assertTrue(self.mesh.dim, 2) + + def test_mesh_nc(self): + self.assertTrue(np.all(self.mesh.n == [6, 2])) + + def test_mesh_nc_xyz(self): + x = np.all(self.mesh.nCx == 6) + y = np.all(self.mesh.nCy == 2) + z = self.mesh.nCz is None + + self.assertTrue(np.all([x, y, z])) + + def test_mesh_nf(self): + x = np.all(self.mesh.nFx == [7, 2]) + y = np.all(self.mesh.nFy == [6, 3]) + z = self.mesh.nFz is None + + self.assertTrue(np.all([x, y, z])) + + def test_mesh_ne(self): + x = np.all(self.mesh.nEx == [6, 3]) + y = np.all(self.mesh.nEy == [7, 2]) + z = self.mesh.nEz is None + + self.assertTrue(np.all([x, y, z])) + + def test_mesh_numbers(self): + c = self.mesh.nC == 12 + f = np.all(self.mesh.nF == [14, 18]) + e = np.all(self.mesh.nE == [18, 14]) + + self.assertTrue(np.all([c, f, e])) + +if __name__ == '__main__': + unittest.main() diff --git a/code/tests/test_mesh.py b/code/tests/test_mesh.py deleted file mode 100644 index a700391f..00000000 --- a/code/tests/test_mesh.py +++ /dev/null @@ -1,74 +0,0 @@ -import unittest -import sys -sys.path.append('../') -from Mesh import Mesh -import numpy as np - - -class TestMeshNumbers3D(unittest.TestCase): - - def setUp(self): - self.mesh = Mesh([6, 2, 3]) - - def test_meshDimensions(self): - self.assertTrue(self.mesh.dim, 3) - - def test_mesh_nc(self): - self.assertTrue(np.all(self.mesh.nc == [6, 2, 3])) - - def test_mesh_nf(self): - x = np.all(self.mesh.nfx == [7, 2, 3]) - y = np.all(self.mesh.nfy == [6, 3, 3]) - z = np.all(self.mesh.nfz == [6, 2, 4]) - - self.assertTrue(np.all([x, y, z])) - - def test_mesh_ne(self): - x = np.all(self.mesh.nex == [6, 3, 4]) - y = np.all(self.mesh.ney == [7, 2, 4]) - z = np.all(self.mesh.nez == [7, 3, 3]) - - self.assertTrue(np.all([x, y, z])) - - def test_mesh_numbers(self): - c = self.mesh.ncells == 36 - f = np.all(self.mesh.nfaces == [42, 54, 48]) - e = np.all(self.mesh.nedges == [72, 56, 63]) - - self.assertTrue(np.all([c, f, e])) - - -class TestMeshNumbers2D(unittest.TestCase): - - def setUp(self): - self.mesh = Mesh([6, 2]) - - def test_meshDimensions(self): - self.assertTrue(self.mesh.dim, 2) - - def test_mesh_nc(self): - self.assertTrue(np.all(self.mesh.nc == [6, 2])) - - def test_mesh_nf(self): - x = np.all(self.mesh.nfx == [7, 2]) - y = np.all(self.mesh.nfy == [6, 3]) - z = self.mesh.nfz is None - - self.assertTrue(np.all([x, y, z])) - - def test_mesh_ne(self): - x = np.all(self.mesh.nex == [6, 3]) - y = np.all(self.mesh.ney == [7, 2]) - z = self.mesh.nez is None - - self.assertTrue(np.all([x, y, z])) - - def test_mesh_numbers(self): - c = self.mesh.ncells == 12 - f = np.all(self.mesh.nfaces == [14, 18]) - e = np.all(self.mesh.nedges == [18, 14]) - - self.assertTrue(np.all([c, f, e])) - -if __name__ == '__main__': - unittest.main() From 6c4910a15733e7735be6e35ce1a17124045721d6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 9 Jul 2013 19:46:50 -0700 Subject: [PATCH 3/6] Base Mesh Documentation --- code/BaseMesh.py | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/code/BaseMesh.py b/code/BaseMesh.py index 90ad07f3..2b72ee80 100644 --- a/code/BaseMesh.py +++ b/code/BaseMesh.py @@ -2,7 +2,33 @@ import numpy as np class BaseMesh(object): - """BaseMesh does all the counting you don't want to do.""" + """BaseMesh does all the counting you don't want to do. + + x0 origin ndarray (dim, ) + n number of cells ndarray (dim, ) + dim dimension of mesh int 1, 2, or 3 + + nCx num cells in x dir int + nCy num cells in y dir int + nCz num cells in z dir int + nC total number of cells int + + nNx num nodes in x dir int + nNy num nodes in y dir int + nNz num nodes in z dir int + nN total number of nodes int + + nEx num edges in x dir ndarray [nEx_x, nEx_y, nEx_z] + nEy num edges in y dir ndarray [nEy_x, nEy_y, nEy_z] + nEz num edges in z dir ndarray [nEz_x, nEz_y, nEz_z] + nE total number of edges ndarray (dim, ) + + nFx num faces in x dir ndarray [nFx_x, nFx_y, nFx_z] + nFy num faces in y dir ndarray [nFy_x, nFy_y, nFy_z] + nFz num faces in z dir ndarray [nFz_x, nFz_y, nFz_z] + nF total number of faces ndarray (dim, ) + + """ def __init__(self, n, x0=None): # Check inputs From 8e3e0e0faae1d2d4072493894f4ae80e05c60cbb Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 9 Jul 2013 19:50:40 -0700 Subject: [PATCH 4/6] TensorMesh now inherits BaseMesh (worked with Luz!) tests for tensorMesh and utils (e.g. ndgrid) are included and pass Split the TensorMesh into Grid and View --- code/TensorGrid.py | 144 ++++++++++++++++ code/TensorMesh.py | 313 ++++++++-------------------------- code/TensorView.py | 96 +++++++++++ code/tests/test_tensorMesh.py | 34 ++++ code/tests/test_utils.py | 53 ++++++ code/utils.py | 64 ++++--- 6 files changed, 439 insertions(+), 265 deletions(-) create mode 100644 code/TensorGrid.py create mode 100644 code/TensorView.py create mode 100644 code/tests/test_tensorMesh.py create mode 100644 code/tests/test_utils.py diff --git a/code/TensorGrid.py b/code/TensorGrid.py new file mode 100644 index 00000000..aa4b31fd --- /dev/null +++ b/code/TensorGrid.py @@ -0,0 +1,144 @@ +import numpy as np +from utils import ndgrid + + +class TensorGrid(object): + """ + Define nodal, cell-centered and staggered tensor grids for 1, 2 and 3 + dimensions. + + This class is inherited by TensorMesh + """ + def __init__(self): + pass + + def vectorNx(): + doc = "Nodal grid vector (1D) in the x direction." + fget = lambda self: np.r_[0., self.hx.cumsum()] + self.x0[0] + return locals() + vectorNx = property(**vectorNx()) + + def vectorNy(): + doc = "Nodal grid vector (1D) in the y direction." + fget = lambda self: None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1] + return locals() + vectorNy = property(**vectorNy()) + + def vectorNz(): + doc = "Nodal grid vector (1D) in the z direction." + fget = lambda self: None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2] + return locals() + vectorNz = property(**vectorNz()) + + def vectorCCx(): + doc = "Cell-centered grid vector (1D) in the x direction." + fget = lambda self: np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0] + return locals() + vectorCCx = property(**vectorCCx()) + + def vectorCCy(): + doc = "Cell-centered grid vector (1D) in the y direction." + fget = lambda self: None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1] + return locals() + vectorCCy = property(**vectorCCy()) + + def vectorCCz(): + doc = "Cell-centered grid vector (1D) in the z direction." + fget = lambda self: None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2] + return locals() + vectorCCz = property(**vectorCCz()) + + def gridCC(): + doc = "Cell-centered grid." + + def fget(self): + if self._gridCC is None: + self._gridCC = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorCCz] if not x is None]) + return self._gridCC + return locals() + _gridCC = None # Store grid by default + gridCC = property(**gridCC()) + + def gridN(): + doc = "Nodal grid." + + def fget(self): + if self._gridN is None: + self._gridN = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorNz] if not x is None]) + return self._gridN + return locals() + _gridN = None # Store grid by default + gridN = property(**gridN()) + + def gridFx(): + doc = "Face staggered grid in the x direction." + + def fget(self): + if self._gridFx is None: + self._gridFx = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorCCz] if not x is None]) + return self._gridFx + return locals() + _gridFx = None # Store grid by default + gridFx = property(**gridFx()) + + def gridFy(): + doc = "Face staggered grid in the y direction." + + def fget(self): + if self._gridFy is None: + self._gridFy = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorCCz] if not x is None]) + return self._gridFy + return locals() + _gridFy = None # Store grid by default + gridFy = property(**gridFy()) + + def gridFz(): + doc = "Face staggered grid in the z direction." + + def fget(self): + if self._gridFz is None: + self._gridFz = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorNz] if not x is None]) + return self._gridFz + return locals() + _gridFz = None # Store grid by default + gridFz = property(**gridFz()) + + def gridEx(): + doc = "Edge staggered grid in the x direction." + + def fget(self): + if self._gridEx is None: + self._gridEx = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorNz] if not x is None]) + return self._gridEx + return locals() + _gridEx = None # Store grid by default + gridEx = property(**gridEx()) + + def gridEy(): + doc = "Edge staggered grid in the y direction." + + def fget(self): + if self._gridEy is None: + self._gridEy = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorNz] if not x is None]) + return self._gridEy + return locals() + _gridEy = None # Store grid by default + gridEy = property(**gridEy()) + + def gridEz(): + doc = "Edge staggered grid in the z direction." + + def fget(self): + if self._gridEz is None: + self._gridEz = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorCCz] if not x is None]) + return self._gridEz + return locals() + _gridEz = None # Store grid by default + gridEz = property(**gridEz()) + + def getBoundaryIndex(self, gridType): + """Needed for faces edges and cells""" + pass + + def getCellNumbering(self): + pass diff --git a/code/TensorMesh.py b/code/TensorMesh.py index 3a061561..6b2305f5 100644 --- a/code/TensorMesh.py +++ b/code/TensorMesh.py @@ -1,257 +1,82 @@ -#------------------------------------------------------------------------------- -# Packages -#------------------------------------------------------------------------------- import numpy as np -import matplotlib.pyplot as plt -from mpl_toolkits.mplot3d import Axes3D +from BaseMesh import BaseMesh +from TensorGrid import TensorGrid +from TensorView import TensorView -#------------------------------------------------------------------------------- -# Class definition -#------------------------------------------------------------------------------- -class TensorMesh: - """ - Define nodal, cell-centered and staggered tensor meshes for 1, 2 and 3 - dimensions. + +class TensorMesh(BaseMesh, TensorGrid, TensorView): """ - - # ---------------------- Properties -------------------------------------- - h = None # Array Spacing or cell-sizes in each direction - x0 = None # Array Origin (x1,x2,x3) - dim = None # Int Dimension - n = None # Array Number of cells in each direction - nC = None # Int Total number of cells - nE = None # Int Total number of edges - nF = None # Int Total number of faces - - # ----------------------- Methods ---------------------------------------- - def __init__(self,h,x0): - """Compute number of edges,faces and cell-centers """ - - # Assign values to properties - self.h = h - self.x0 = x0 - - - #Compute derived properties - self.dim = np.size(x0) - - #Compute the num of cells in each direction - self.n = np.zeros((self.dim,1)) - - for d in range(self.dim): - self.n[d] = np.size(h[d]) - - # Compute the number of cell-centers - self.nC = np.prod(self.n) - - # Compute the number of edges (makes sense only for 3D) - # Equivalent to: - # nEdges = n[0] * (n[1]+1) * (n[2]+1) - # + (n[0]+1) * ny[1] * (nz[2]+1) - # + (n[0]+1) * (ny[1]+1)* (nz[2]) - - if self.dim == 3: - self.nE = np.prod(np.kron(np.ones((3,1)),self.n.T)+np.ones((3,3))-np.eye(3),1) - print self.nE - - # Compute the number of faces (makes sense only for 2 and 3D) - # Equivalent to - # nFaces = (n[0]+1) * n[1] * n[2] - # + n[0] * (ny[1]+1) * nz[2] - # + n[0] * ny[1] * (nz[2]+1) - if self.dim >=2: - self.nF = np.prod(np.kron(np.ones((self.dim,1)),self.n.T)+np.eye(self.dim),1) - print self.nF + TensorMesh is a mesh class that deals with tensor product meshes. - def xin(self,i): - """Construct the 1D nodal mesh from the ith-component of h. Return an array.""" - return np.insert(np.cumsum(self.h[i-1]),0,0.0) + self.x0[i-1] - - def xic(self,i): - """Construct the 1D cell-centerd mesh from the ith-component of h. Return an array.""" - return .5*( np.insert(np.cumsum(self.h[i-1][:,0:-1]),0,0.0) + np.cumsum(self.h[i-1])) + Any Mesh that has a constant width along the entire axis + such that it can defined by a single width vector, called 'h'. + + e.g. + + hx = np.array([1,1,1]) + hy = np.array([1,2]) + hz = np.array([1,1,1,1]) + + mesh = TensorMesh([hx, hy, hz]) + + """ + def __init__(self, h, x0=None): + super(TensorMesh, self).__init__(np.array([len(x) for x in h]), x0) + + assert len(h) == len(x0), "Dimension mismatch. x0 != len(h)" + + for i, h_i in enumerate(h): + assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) + + # Ensure h contains 1D vectors + self._h = [x.ravel() for x in h] + + def h(): + doc = "h is a list containing the cell widths of the tensor mesh in each dimension." + fget = lambda self: self._h + return locals() + h = property(**h()) + + def hx(): + doc = "Width of cells in the x direction" + fget = lambda self: self._h[0] + return locals() + hx = property(**hx()) + + def hy(): + doc = "Width of cells in the y direction" + fget = lambda self: None if self.dim < 2 else self._h[1] + return locals() + hy = property(**hy()) + + def hz(): + doc = "Width of cells in the z direction" + fget = lambda self: None if self.dim < 3 else self._h[2] + return locals() + hz = property(**hz()) - def getNodalGrid(self): - """Construct nodal grid for 1, 2 and 3 dimensions""" - if self.dim==1: - return [self.xin(1)] - elif self.dim==2: - return self.ndgrid([self.xin(1),self.xin(2)]) - elif self.dim==3: - return self.ndgrid([self.xin(1),self.xin(2),self.xin(3)]) - - def ndgrid(self, xin): - """Form tensorial grid for 1, 2 and 3 dimensions. Return X1,X2,X3 arrays depending on the dimension""" - ei = lambda i : np.ones((np.size(xin[i-1]),1)) - - if self.dim==1: - return [xin] - elif self.dim==2: - X1 = np.kron(ei(2),xin[0]).reshape(-1,1) - X2 = np.kron(xin[1],ei(1).T).reshape(-1,1) - return X1,X2 - elif self.dim==3: - X1 = np.kron(ei(3),np.kron(ei(2),xin[0])).reshape(-1,1) - X2 = np.kron(ei(3).T,np.kron(xin[1],ei(1).T)).reshape(-1,1) - X3 = np.kron(xin[2],np.kron(ei(2),ei(1))).T.reshape(-1,1) - - return X1,X2,X3 - - def getCellCenteredGrid(self): - """Construct cell-centered grid for 1, 2 and 3 dimensions.""" - - if self.dim==1: - return [self.xic(1)] - elif self.dim==2: - return self.ndgrid([self.xic(1),self.xic(2)]) - elif self.dim==3: - return self.ndgrid([self.xic(1),self.xic(2),self.xic(3)]) - - def getFaceStgGrid(self,direction): - """Construct the face staggered grids for 2 and 3 dimensions.""" - if self.dim==1: - print 'Error: dimension must be larger than 1' - elif self.dim==2: - if direction == 1: - return self.ndgrid([self.xin(1),self.xic(2)]) - elif direction == 2: - return self.ndgrid([self.xic(1),self.xin(2)]) - else: - print 'Error: direction must be equal to 1 or 2' - elif self.dim==3: - if direction == 1: - return self.ndgrid([self.xin(1),self.xic(2),self.xic(3)]) - elif direction == 2: - return self.ndgrid([self.xic(1),self.xin(2),self.xic(3)]) - elif direction == 3: - return self.ndgrid([self.xic(1),self.xic(2),self.xin(3)]) - else: - print 'Error: direction must be equal to 1, 2 or 3' - - def getEdgeStgGrid(self,direction): - """Construct the edge staggered grids for 3 dimension case.""" - if self.dim != 3: - print 'Error: dimension must be equal to 3' - else: - if direction == 1: - return self.ndgrid([self.xic(1),self.xin(2),self.xin(3)]) - elif direction == 2: - return self.ndgrid([self.xin(1),self.xic(2),self.xin(3)]) - elif direction == 3: - return self.ndgrid([self.xin(1),self.xin(2),self.xic(3)]) - else: - print 'Error: direction must be equal to 1, 2 or 3' - - def plotImage(self,I): - - if self.dim==1: - fig = plt.figure(1) - fig.clf() - ax=plt.subplot(111) - if np.size(I)==self.n[0]: - print 'cell-centered image' - xx = self.getCellCenteredGrid() - ax.plot(xx[0],I,'ro') - elif np.size(I)==self.n[0]+1: - print 'nodal image' - xx = self.getNodalGrid() - ax.plot(xx[0],I,'bs') - - fig.show() - - - def plotGrid(self): - """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.""" - if self.dim == 1: - fig = plt.figure(1) - fig.clf() - ax = plt.subplot(111) - xn = self.getNodalGrid() - xc = self.getCellCenteredGrid() - print xn - ax.hold(True) - ax.plot(xn,np.ones(np.shape(xn)),'bs') - ax.plot(xc,np.ones(np.shape(xc)),'ro') - ax.plot(xn,np.ones(np.shape(xn)),'k--') - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - fig.show() - elif self.dim == 2: - fig = plt.figure(2) - fig.clf() - ax = plt.subplot(111) - xn = self.getNodalGrid() - xc = self.getCellCenteredGrid() - xs1 = self.getFaceStgGrid(1) - xs2 = self.getFaceStgGrid(2) - - ax.hold(True) - ax.plot(xn[0],xn[1],'bs') - ax.plot(xc[0],xc[1],'ro') - ax.plot(xs1[0],xs1[1],'g>') - ax.plot(xs2[0],xs2[1],'g^') - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - fig.show() - elif self.dim == 3: - fig = plt.figure(3) - fig.clf() - ax = fig.add_subplot(111, projection='3d') - xn = self.getNodalGrid() - xc = self.getCellCenteredGrid() - xfs1 = self.getFaceStgGrid(1) - xfs2 = self.getFaceStgGrid(2) - xfs3 = self.getFaceStgGrid(3) - xes1 = self.getEdgeStgGrid(1) - xes2 = self.getEdgeStgGrid(2) - xes3 = self.getEdgeStgGrid(3) - - ax.hold(True) - ax.plot(xn[0],xn[1],'bs',zs=xn[2]) - ax.plot(xc[0],xc[1],'ro',zs=xc[2]) - ax.plot(xfs1[0],xfs1[1],'g>',zs=xfs1[2]) - ax.plot(xfs2[0],xfs2[1],'g<',zs=xfs2[2]) - ax.plot(xfs3[0],xfs3[1],'g^',zs=xfs3[2]) - ax.plot(xes1[0],xes1[1],'k>',zs=xes1[2]) - ax.plot(xes2[0],xes2[1],'k<',zs=xes2[2]) - ax.plot(xes3[0],xes3[1],'k^',zs=xes3[2]) - ax.grid(True) - ax.hold(False) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - ax.set_zlabel('x3') - fig.show() - - - if __name__ == '__main__': print('Welcome to tensor mesh!') - + testDim = 1 - h1 = 0.3*np.ones((1,7)) - h1[:,0] = 0.5 - h1[:,-1] = 0.6 - h2 = .5* np.ones((1,4)) - h3 = .4* np.ones((1,6)) - x0 = np.zeros((3,1)) - + h1 = 0.3*np.ones((1, 7)) + h1[:, 0] = 0.5 + h1[:, -1] = 0.6 + h2 = .5 * np.ones((1, 4)) + h3 = .4 * np.ones((1, 6)) + x0 = np.zeros((3, 1)) + if testDim == 1: h = [h1] - x0 = x0[0] - elif testDim==2: - h = [h1,h2] + x0 = x0[0] + elif testDim == 2: + h = [h1, h2] x0 = x0[0:2] else: - h = [h1,h2,h3] - - I = np.linspace(0,1,8) - M = TensorMesh(h,x0) - - xn = M.plotGrid() - - - \ No newline at end of file + h = [h1, h2, h3] + + I = np.linspace(0, 1, 8) + M = TensorMesh(h, x0) + + xn = M.plotGrid() diff --git a/code/TensorView.py b/code/TensorView.py new file mode 100644 index 00000000..18a7bd00 --- /dev/null +++ b/code/TensorView.py @@ -0,0 +1,96 @@ +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + + +class TensorView(object): + """ + Provides viewing functions for TensorMesh + + This class is inherited by TensorMesh + """ + def __init__(self): + pass + + def plotImage(self, I): + + if self.dim == 1: + fig = plt.figure(1) + fig.clf() + ax = plt.subplot(111) + if np.size(I) == self.n[0]: + print 'cell-centered image' + xx = self.gridCC + ax.plot(xx[0], I, 'ro') + elif np.size(I) == self.n[0]+1: + print 'nodal image' + xx = self.gridN + ax.plot(xx[0], I, 'bs') + + fig.show() + + def plotGrid(self): + """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.""" + if self.dim == 1: + fig = plt.figure(1) + fig.clf() + ax = plt.subplot(111) + xn = self.gridN + xc = self.gridCC + print xn + ax.hold(True) + ax.plot(xn, np.ones(np.shape(xn)), 'bs') + ax.plot(xc, np.ones(np.shape(xc)), 'ro') + ax.plot(xn, np.ones(np.shape(xn)), 'k--') + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + fig.show() + elif self.dim == 2: + fig = plt.figure(2) + fig.clf() + ax = plt.subplot(111) + xn = self.gridN + xc = self.gridCC + xs1 = self.gridFx + xs2 = self.gridFy + + ax.hold(True) + ax.plot(xn[:, 0], xn[:, 1], 'bs') + ax.plot(xc[:, 0], xc[:, 1], 'ro') + ax.plot(xs1[:, 0], xs1[:, 1], 'g>') + ax.plot(xs2[:, 0], xs2[:, 1], 'g^') + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + fig.show() + elif self.dim == 3: + fig = plt.figure(3) + fig.clf() + ax = fig.add_subplot(111, projection='3d') + xn = self.gridN + xc = self.gridCC + xfs1 = self.gridFx + xfs2 = self.gridFy + xfs3 = self.gridFz + + xes1 = self.gridEx + xes2 = self.gridEy + xes3 = self.gridEz + + ax.hold(True) + ax.plot(xn[:, 0], xn[:, 1], 'bs', zs=xn[:, 2]) + ax.plot(xc[:, 0], xc[:, 1], 'ro', zs=xc[:, 2]) + ax.plot(xfs1[:, 0], xfs1[:, 1], 'g>', zs=xfs1[:, 2]) + ax.plot(xfs2[:, 0], xfs2[:, 1], 'g<', zs=xfs2[:, 2]) + ax.plot(xfs3[:, 0], xfs3[:, 1], 'g^', zs=xfs3[:, 2]) + ax.plot(xes1[:, 0], xes1[:, 1], 'k>', zs=xes1[:, 2]) + ax.plot(xes2[:, 0], xes2[:, 1], 'k<', zs=xes2[:, 2]) + ax.plot(xes3[:, 0], xes3[:, 1], 'k^', zs=xes3[:, 2]) + ax.grid(True) + ax.hold(False) + ax.set_xlabel('x1') + ax.set_ylabel('x2') + ax.set_zlabel('x3') + fig.show() diff --git a/code/tests/test_tensorMesh.py b/code/tests/test_tensorMesh.py new file mode 100644 index 00000000..e1b74374 --- /dev/null +++ b/code/tests/test_tensorMesh.py @@ -0,0 +1,34 @@ +import numpy as np +import unittest +import sys +sys.path.append('../') +from TensorMesh import TensorMesh + + +class TestSequenceFunctions(unittest.TestCase): + + def setUp(self): + a = np.array([1, 1, 1]) + b = np.array([1, 2]) + x0 = np.array([3, 5]) + self.mesh2 = TensorMesh([a, b], x0) + + 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) + + def test_vectorCC_2D(self): + testNx = np.array([3.5, 4.5, 5.5]) + testNy = np.array([5.5, 7]) + + xtest = np.all(self.mesh2.vectorCCx == testNx) + ytest = np.all(self.mesh2.vectorCCy == testNy) + self.assertTrue(xtest and ytest) + + +if __name__ == '__main__': + unittest.main() diff --git a/code/tests/test_utils.py b/code/tests/test_utils.py new file mode 100644 index 00000000..5d05b710 --- /dev/null +++ b/code/tests/test_utils.py @@ -0,0 +1,53 @@ +import numpy as np +import unittest +import sys +sys.path.append('../') +from utils import mkvc, ndgrid + + +class TestSequenceFunctions(unittest.TestCase): + + def setUp(self): + self.a = np.array([1, 2, 3]) + self.b = np.array([1, 2]) + self.c = np.array([1, 2, 3, 4]) + + def test_mkvc1(self): + x = mkvc(self.a) + self.assertTrue(x.shape, (3,)) + + def test_mkvc2(self): + x = mkvc(self.a, 2) + self.assertTrue(x.shape, (3, 1)) + + def test_mkvc3(self): + x = mkvc(self.a, 3) + self.assertTrue(x.shape, (3, 1, 1)) + + def test_ndgrid_2D(self): + XY = ndgrid([self.a, self.b]) + + X1_test = np.array([1, 2, 3, 1, 2, 3]) + X2_test = np.array([1, 1, 1, 2, 2, 2]) + + xtest = np.all(XY[:, 0] == X1_test) + ytest = np.all(XY[:, 1] == X2_test) + + self.assertTrue(xtest and ytest) + + def test_ndgrid_3D(self): + XYZ = ndgrid([self.a, self.b, self.c]) + + X1_test = np.array([1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3]) + X2_test = np.array([1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) + X3_test = np.array([1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4]) + + xtest = np.all(XYZ[:, 0] == X1_test) + ytest = np.all(XYZ[:, 1] == X2_test) + ztest = np.all(XYZ[:, 2] == X3_test) + + self.assertTrue(xtest and ytest and ztest) + +if __name__ == '__main__': + unittest.main() + diff --git a/code/utils.py b/code/utils.py index 7ea54f17..a41d8687 100644 --- a/code/utils.py +++ b/code/utils.py @@ -1,4 +1,5 @@ from numpy import * +import numpy as np def diff(A, d): @@ -43,35 +44,56 @@ def ave(A, d): print('d must be 1,2 or 3') -def reshapeF(sp, d): - return reshape(sp, d, 'F') +def reshapeF(x, size): + return np.reshape(x, size, order='F') -def mkvc(A): - return reshape(A, [size(A), 1], 'F').flatten() +def mkvc(x, numDims=1): + """Creates a vector with the number of dimension specified + + e.g.: + + a = np.array(1,2,3) + + mkvc(a, 1).shape + > (3, ) + + mkvc(a, 2).shape + > (3, 1) + + mkvc(a, 3).shape + > (3, 1, 1) + + """ + assert type(x) == np.ndarray, "Vector must be a numpy array" + + if numDims == 1: + return x.flatten(order='F') + elif numDims == 2: + return x.flatten(order='F')[:, np.newaxis] + elif numDims == 3: + return x.flatten(order='F')[:, np.newaxis, np.newaxis] -def ndgrid(x, y, z): +def ndgrid(xin): + """Form tensorial grid for 1, 2 and 3 dimensions. Return X1,X2,X3 arrays depending on the dimension""" - n1 = size(x) - n2 = size(y) - n3 = size(z) - X = zeros([n1, n2, n3]) - Y = zeros([n1, n2, n3]) - Z = zeros([n1, n2, n3]) - for i in range(0, n2): - for j in range(0, n3): - X[:, i, j] = x + if len(xin) == 1: + return xin + elif len(xin) == 2: + X2, X1 = [mkvc(x) for x in np.broadcast_arrays(mkvc(xin[1], 1), mkvc(xin[0], 2))] + return np.c_[X1, X2] + elif len(xin) == 3: + X3, X2, X1 = [mkvc(x) for x in np.broadcast_arrays(mkvc(xin[2], 1), mkvc(xin[1], 2), mkvc(xin[0], 3))] + return np.c_[X1, X2, X3] - for i in range(0, n1): - for j in range(0, n3): - Y[i, :, j] = y - for i in range(0, n1): - for j in range(0, n2): - Z[i, j, :] = z +def flattenF(x): + return np.flatten(x, order='F') - return (X, Y, Z) + +def printF(x): + pass def ind2sub(shape, ind): From e2fa98602c1bde796c1606279013978eb32c2610 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Tue, 9 Jul 2013 21:38:17 -0700 Subject: [PATCH 5/6] Moved code not being used in new code to it's own directory. --- code/{ => EldadsCode}/GaussNewton.py | 0 code/{ => EldadsCode}/MFDdriver.py | 0 code/{ => EldadsCode}/getCellVolume.py | 0 code/{ => EldadsCode}/getDiffOps.py | 0 code/{ => EldadsCode}/getEdgeInnerProduct.py | 0 code/{ => EldadsCode}/getEdgeTangent.py | 0 code/{ => EldadsCode}/getFaceInnerProduct.py | 0 code/{ => EldadsCode}/getFaceNormals.py | 0 code/{ => EldadsCode}/getVolume.py | 0 code/{ => EldadsCode}/inv3X3BlockDiagonal.py | 0 code/{ => EldadsCode}/meshUtils.py | 0 code/{ => EldadsCode}/sputils.py | 0 code/{ => EldadsCode}/tools.py | 0 code/{ => EldadsCode}/zevel.py | 0 14 files changed, 0 insertions(+), 0 deletions(-) rename code/{ => EldadsCode}/GaussNewton.py (100%) rename code/{ => EldadsCode}/MFDdriver.py (100%) rename code/{ => EldadsCode}/getCellVolume.py (100%) rename code/{ => EldadsCode}/getDiffOps.py (100%) rename code/{ => EldadsCode}/getEdgeInnerProduct.py (100%) rename code/{ => EldadsCode}/getEdgeTangent.py (100%) rename code/{ => EldadsCode}/getFaceInnerProduct.py (100%) rename code/{ => EldadsCode}/getFaceNormals.py (100%) rename code/{ => EldadsCode}/getVolume.py (100%) rename code/{ => EldadsCode}/inv3X3BlockDiagonal.py (100%) rename code/{ => EldadsCode}/meshUtils.py (100%) rename code/{ => EldadsCode}/sputils.py (100%) rename code/{ => EldadsCode}/tools.py (100%) rename code/{ => EldadsCode}/zevel.py (100%) diff --git a/code/GaussNewton.py b/code/EldadsCode/GaussNewton.py similarity index 100% rename from code/GaussNewton.py rename to code/EldadsCode/GaussNewton.py diff --git a/code/MFDdriver.py b/code/EldadsCode/MFDdriver.py similarity index 100% rename from code/MFDdriver.py rename to code/EldadsCode/MFDdriver.py diff --git a/code/getCellVolume.py b/code/EldadsCode/getCellVolume.py similarity index 100% rename from code/getCellVolume.py rename to code/EldadsCode/getCellVolume.py diff --git a/code/getDiffOps.py b/code/EldadsCode/getDiffOps.py similarity index 100% rename from code/getDiffOps.py rename to code/EldadsCode/getDiffOps.py diff --git a/code/getEdgeInnerProduct.py b/code/EldadsCode/getEdgeInnerProduct.py similarity index 100% rename from code/getEdgeInnerProduct.py rename to code/EldadsCode/getEdgeInnerProduct.py diff --git a/code/getEdgeTangent.py b/code/EldadsCode/getEdgeTangent.py similarity index 100% rename from code/getEdgeTangent.py rename to code/EldadsCode/getEdgeTangent.py diff --git a/code/getFaceInnerProduct.py b/code/EldadsCode/getFaceInnerProduct.py similarity index 100% rename from code/getFaceInnerProduct.py rename to code/EldadsCode/getFaceInnerProduct.py diff --git a/code/getFaceNormals.py b/code/EldadsCode/getFaceNormals.py similarity index 100% rename from code/getFaceNormals.py rename to code/EldadsCode/getFaceNormals.py diff --git a/code/getVolume.py b/code/EldadsCode/getVolume.py similarity index 100% rename from code/getVolume.py rename to code/EldadsCode/getVolume.py diff --git a/code/inv3X3BlockDiagonal.py b/code/EldadsCode/inv3X3BlockDiagonal.py similarity index 100% rename from code/inv3X3BlockDiagonal.py rename to code/EldadsCode/inv3X3BlockDiagonal.py diff --git a/code/meshUtils.py b/code/EldadsCode/meshUtils.py similarity index 100% rename from code/meshUtils.py rename to code/EldadsCode/meshUtils.py diff --git a/code/sputils.py b/code/EldadsCode/sputils.py similarity index 100% rename from code/sputils.py rename to code/EldadsCode/sputils.py diff --git a/code/tools.py b/code/EldadsCode/tools.py similarity index 100% rename from code/tools.py rename to code/EldadsCode/tools.py diff --git a/code/zevel.py b/code/EldadsCode/zevel.py similarity index 100% rename from code/zevel.py rename to code/EldadsCode/zevel.py From 68b4719aea6092e9ab8f174f7a11c246a6a60aa6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 9 Jul 2013 22:00:17 -0700 Subject: [PATCH 6/6] Ability to run all tests at once. Removed __init__.py added runTests.py and .sh so that we can run all tests at one time. --- code/tests/__init__.py | 3 --- code/tests/runTests.py | 11 +++++++++++ code/tests/runTests.sh | 3 +++ code/tests/test_basemesh.py | 2 +- code/utils.py | 16 +++++++--------- 5 files changed, 22 insertions(+), 13 deletions(-) delete mode 100644 code/tests/__init__.py create mode 100644 code/tests/runTests.py create mode 100644 code/tests/runTests.sh diff --git a/code/tests/__init__.py b/code/tests/__init__.py deleted file mode 100644 index 80b21e43..00000000 --- a/code/tests/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -# init.py - -from test_mesh import * \ No newline at end of file diff --git a/code/tests/runTests.py b/code/tests/runTests.py new file mode 100644 index 00000000..b6662f00 --- /dev/null +++ b/code/tests/runTests.py @@ -0,0 +1,11 @@ +import glob +import unittest + +# This code will run all tests in directory named test_*.py + +test_file_strings = glob.glob('test_*.py') +module_strings = [str[0:len(str)-3] for str in test_file_strings] +suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] +testSuite = unittest.TestSuite(suites) +text_runner = unittest.TextTestRunner().run(testSuite) diff --git a/code/tests/runTests.sh b/code/tests/runTests.sh new file mode 100644 index 00000000..9bd83e67 --- /dev/null +++ b/code/tests/runTests.sh @@ -0,0 +1,3 @@ +#!/bin/sh + +python -m unittest discover \ No newline at end of file diff --git a/code/tests/test_basemesh.py b/code/tests/test_basemesh.py index 88655c60..9aece19b 100644 --- a/code/tests/test_basemesh.py +++ b/code/tests/test_basemesh.py @@ -5,7 +5,7 @@ from BaseMesh import BaseMesh import numpy as np -class TestMeshNumbers3D(unittest.TestCase): +class TestBaseMesh(unittest.TestCase): def setUp(self): self.mesh = BaseMesh([6, 2, 3]) diff --git a/code/utils.py b/code/utils.py index a41d8687..a5b877cb 100644 --- a/code/utils.py +++ b/code/utils.py @@ -75,9 +75,15 @@ def mkvc(x, numDims=1): return x.flatten(order='F')[:, np.newaxis, np.newaxis] -def ndgrid(xin): +def ndgrid(*args): """Form tensorial grid for 1, 2 and 3 dimensions. Return X1,X2,X3 arrays depending on the dimension""" + # you can either pass a list [x1, x2, x3] or each seperately + if type(args[0]) == list: + xin = args[0] + else: + xin = args + if len(xin) == 1: return xin elif len(xin) == 2: @@ -88,14 +94,6 @@ def ndgrid(xin): return np.c_[X1, X2, X3] -def flattenF(x): - return np.flatten(x, order='F') - - -def printF(x): - pass - - def ind2sub(shape, ind): # From the given shape, returns the subscrips of the given index revshp = []