From 3ad2c98d43baf42bbaf94ef1aadf6103c43ecb70 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 21 Jul 2016 09:58:14 -0700 Subject: [PATCH] use @property decorator in curve mesh (see: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d?groups=code_patterns%3A3bECxdfc%3Af0&tab=basics) --- SimPEG/Mesh/CurvilinearMesh.py | 435 ++++++++++++++++------------- tests/mesh/test_CurvilinearMesh.py | 96 +++++-- 2 files changed, 300 insertions(+), 231 deletions(-) diff --git a/SimPEG/Mesh/CurvilinearMesh.py b/SimPEG/Mesh/CurvilinearMesh.py index abc9ff3a..d3c0c5b9 100644 --- a/SimPEG/Mesh/CurvilinearMesh.py +++ b/SimPEG/Mesh/CurvilinearMesh.py @@ -4,14 +4,28 @@ from DiffOperators import DiffOperators from InnerProducts import InnerProducts from View import CurvView + # Some helper functions. -length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 -length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5 -normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2)) -normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2)) +def length2D(x): + return (x[:, 0]**2 + x[:, 1]**2)**0.5 -class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView): +def length3D(x): + return (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5 + + +def normalize2D(x): + return x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2)) + + +def normalize3D(x): + return x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2)) + + +# Curvi Mesh + +class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, + CurvView): """ CurvilinearMesh is a mesh class that deals with curvilinear meshes. @@ -53,7 +67,10 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie def fget(self): if self._gridCC 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 return locals() _gridCC = None # Store grid by default @@ -70,99 +87,93 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie _gridN = None # Store grid by default gridN = property(**gridN()) - def gridFx(): - doc = "Face staggered grid in the x direction." + @property + def gridFx(self): + """ + Face staggered grid in the x direction. + """ - def fget(self): - if self._gridFx is None: - N = self.r(self.gridN, 'N', 'N', 'M') - if self.dim == 2: - XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] - self._gridFx = np.c_[XY[0], XY[1]] - elif self.dim == 3: - XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N] - self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]] - return self._gridFx - return locals() - _gridFx = None # Store grid by default - gridFx = property(**gridFx()) + if getattr(self, '_gridFx', None) is None: + N = self.r(self.gridN, 'N', 'N', 'M') + if self.dim == 2: + XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] + self._gridFx = np.c_[XY[0], XY[1]] + elif self.dim == 3: + XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N] + self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]] + return self._gridFx - def gridFy(): - doc = "Face staggered grid in the y direction." + @property + def gridFy(self): + """ + Face staggered grid in the y direction. + """ - def fget(self): - if self._gridFy is None: - N = self.r(self.gridN, 'N', 'N', 'M') - if self.dim == 2: - XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] - self._gridFy = np.c_[XY[0], XY[1]] - elif self.dim == 3: - XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N] - self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]] - return self._gridFy - return locals() - _gridFy = None # Store grid by default - gridFy = property(**gridFy()) + if getattr(self, '_gridFy', None) is None: + N = self.r(self.gridN, 'N', 'N', 'M') + if self.dim == 2: + XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] + self._gridFy = np.c_[XY[0], XY[1]] + elif self.dim == 3: + XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + + n[1:, :, :-1] + n[1:, :, 1:])) for n in N] + self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]] + return self._gridFy - def gridFz(): - doc = "Face staggered grid in the z direction." + @property + def gridFz(self): + """ + Face staggered grid in the y direction. + """ - def fget(self): - if self._gridFz is None and self.dim == 3: - N = self.r(self.gridN, 'N', 'N', 'M') - XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N] - self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]] - return self._gridFz - return locals() - _gridFz = None # Store grid by default - gridFz = property(**gridFz()) + if getattr(self, '_gridFz', None) is None: + N = self.r(self.gridN, 'N', 'N', 'M') + XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + + n[1:, :-1, :] + n[1:, 1:, :])) for n in N] + self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]] + return self._gridFz - def gridEx(): - doc = "Edge staggered grid in the x direction." + @property + def gridEx(self): + """ + Edge staggered grid in the x direction. + """ + if getattr(self, '_gridEx', None) is None: + N = self.r(self.gridN, 'N', 'N', 'M') + if self.dim == 2: + XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] + self._gridEx = np.c_[XY[0], XY[1]] + elif self.dim == 3: + XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N] + self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]] + return self._gridEx - def fget(self): - if self._gridEx is None: - N = self.r(self.gridN, 'N', 'N', 'M') - if self.dim == 2: - XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] - self._gridEx = np.c_[XY[0], XY[1]] - elif self.dim == 3: - XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N] - self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]] - return self._gridEx - return locals() - _gridEx = None # Store grid by default - gridEx = property(**gridEx()) + @property + def gridEy(self): + """ + Edge staggered grid in the y direction. + """ + if getattr(self, '_gridEy', None) is None: + N = self.r(self.gridN, 'N', 'N', 'M') + if self.dim == 2: + XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] + self._gridEy = np.c_[XY[0], XY[1]] + elif self.dim == 3: + XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N] + self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]] + return self._gridEy - def gridEy(): - doc = "Edge staggered grid in the y direction." - - def fget(self): - if self._gridEy is None: - N = self.r(self.gridN, 'N', 'N', 'M') - if self.dim == 2: - XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] - self._gridEy = np.c_[XY[0], XY[1]] - elif self.dim == 3: - XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N] - self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]] - 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 and self.dim == 3: - N = self.r(self.gridN, 'N', 'N', 'M') - XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N] - self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]] - return self._gridEz - return locals() - _gridEz = None # Store grid by default - gridEz = property(**gridEz()) + @property + def gridEz(self): + """ + Edge staggered grid in the z direction. + """ + if getattr(self, '_gridEz', None) is None and self.dim == 3: + N = self.r(self.gridN, 'N', 'N', 'M') + XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N] + self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]] + return self._gridEz # --------------- Geometries --------------------- # @@ -194,78 +205,94 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie # | / | / # D -------------- C # node(i+1,j,k) node(i+1,j+1,k) - def vol(): - doc = "Construct cell volumes of the 3D model as 1d array." - def fget(self): - if(self._vol is None): - if self.dim == 2: - A, B, C, D = Utils.indexCube('ABCD', self.vnC+1) - normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D) - self._vol = area - elif self.dim == 3: - # Each polyhedron can be decomposed into 5 tetrahedrons - # However, this presents a choice so we may as well divide in two ways and average. - A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC+1) + @property + def vol(self): + """ + Construct cell volumes of the 3D model as 1d array + """ - vol1 = (Utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top - Utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top - Utils.volTetra(self.gridN, B, D, E, G) + # middle - Utils.volTetra(self.gridN, B, C, D, G) + # cutted edge bottom - Utils.volTetra(self.gridN, D, E, G, H)) # cutted edge bottom + if getattr(self, '_vol', None) is None: + if self.dim == 2: + A, B, C, D = Utils.indexCube('ABCD', self.vnC+1) + normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros( + (self.nN, 1))], A, B, C, D) + self._vol = area + elif self.dim == 3: + # Each polyhedron can be decomposed into 5 tetrahedrons + # However, this presents a choice so we may as well divide in + # two ways and average. + A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC + + 1) - vol2 = (Utils.volTetra(self.gridN, A, F, B, C) + # cutted edge top - Utils.volTetra(self.gridN, A, E, F, H) + # cutted edge top - Utils.volTetra(self.gridN, A, H, F, C) + # middle - Utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom - Utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom + vol1 = (Utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top + Utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top + Utils.volTetra(self.gridN, B, D, E, G) + # middle + Utils.volTetra(self.gridN, B, C, D, G) + # cutted edge bottom + Utils.volTetra(self.gridN, D, E, G, H)) # cutted edge bottom - self._vol = (vol1 + vol2)/2 - return self._vol - return locals() - _vol = None - vol = property(**vol()) + vol2 = (Utils.volTetra(self.gridN, A, F, B, C) + # cutted edge top + Utils.volTetra(self.gridN, A, E, F, H) + # cutted edge top + Utils.volTetra(self.gridN, A, H, F, C) + # middle + Utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom + Utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom - def area(): - doc = "Face areas." + self._vol = (vol1 + vol2)/2 + return self._vol - def fget(self): - if(self._area is None or self._normals is None): - # Compute areas of cell faces - if(self.dim == 2): - xy = self.gridN - A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy])) - edge1 = xy[B, :] - xy[A, :] - normal1 = np.c_[edge1[:, 1], -edge1[:, 0]] - area1 = length2D(edge1) - A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy])) - # Note that we are doing A-D to make sure the normal points the right way. - # Think about it. Look at the picture. Normal points towards C iff you do this. - edge2 = xy[A, :] - xy[D, :] - normal2 = np.c_[edge2[:, 1], -edge2[:, 0]] - area2 = length2D(edge2) - self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] - self._normals = [normalize2D(normal1), normalize2D(normal2)] - elif(self.dim == 3): + @property + def area(self): + if (getattr(self, '_area', None) is None or + getattr(self, '_normals', None) is None): + # Compute areas of cell faces + if(self.dim == 2): + xy = self.gridN + A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, + self.nCy])) + edge1 = xy[B, :] - xy[A, :] + normal1 = np.c_[edge1[:, 1], -edge1[:, 0]] + area1 = length2D(edge1) + A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, + self.nNy])) + # Note that we are doing A-D to make sure the normal points the + # right way. + # Think about it. Look at the picture. Normal points towards C + # iff you do this. + edge2 = xy[A, :] - xy[D, :] + normal2 = np.c_[edge2[:, 1], -edge2[:, 0]] + area2 = length2D(edge2) + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] + self._normals = [normalize2D(normal1), normalize2D(normal2)] - A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array([self.nNx, self.nCy, self.nCz])) - normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False) + elif(self.dim == 3): - A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array([self.nCx, self.nNy, self.nCz])) - normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False) + A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array( + [self.nNx, self.nCy, self.nCz])) + normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, + average=False, + normalizeNormals=False) - A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array([self.nCx, self.nCy, self.nNz])) - normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False) + A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array( + [self.nCx, self.nNy, self.nCz])) + normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, + average=False, + normalizeNormals=False) - self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)] - self._normals = [normal1, normal2, normal3] - return self._area - return locals() - _area = None - area = property(**area()) + A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array( + [self.nCx, self.nCy, self.nNz])) + normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, + average=False, + normalizeNormals=False) - def normals(): - doc = """Face normals: calling this will average + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), + Utils.mkvc(area3)] + self._normals = [normal1, normal2, normal3] + return self._area + + @property + def normals(self): + """ + Face normals: calling this will average the computed normals so that there is one per face. This is especially relevant in 3D, as there are up to 4 different normals @@ -276,58 +303,64 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M') """ - def fget(self): - if(self._normals is None): - self.area # calling .area will create the face normals - if self.dim == 2: - return normalize2D(np.r_[self._normals[0], self._normals[1]]) - elif self.dim == 3: - normal1 = (self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3])/4 - normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4 - normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4 - return normalize3D(np.r_[normal1, normal2, normal3]) - return locals() - _normals = None - normals = property(**normals()) + if getattr(self, '_normals', None) is None: + self.area # calling .area will create the face normals + if self.dim == 2: + return normalize2D(np.r_[self._normals[0], self._normals[1]]) + elif self.dim == 3: + normal1 = (self._normals[0][0] + self._normals[0][1] + self._normals[0][2] + self._normals[0][3])/4 + normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4 + normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4 + return normalize3D(np.r_[normal1, normal2, normal3]) - def edge(): - doc = "Edge legnths." - - def fget(self): - if(self._edge is None or self._tangents is None): - if(self.dim == 2): - xy = self.gridN - A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy])) - edge1 = xy[D, :] - xy[A, :] - A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy])) - edge2 = xy[B, :] - xy[A, :] - self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))] - self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge] - elif(self.dim == 3): - xyz = self.gridN - A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy, self.nNz])) - edge1 = xyz[D, :] - xyz[A, :] - A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy, self.nNz])) - edge2 = xyz[B, :] - xyz[A, :] - A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, self.nNy, self.nCz])) - edge3 = xyz[E, :] - xyz[A, :] - self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))] - self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge] + @property + def edge(self): + """ + Edge lengths + """ + if getattr(self, '_edge', None) is None: + if(self.dim == 2): + xy = self.gridN + A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, + self.nNy])) + edge1 = xy[D, :] - xy[A, :] + A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, + self.nCy])) + edge2 = xy[B, :] - xy[A, :] + self._edge = np.r_[Utils.mkvc(length2D(edge1)), + Utils.mkvc(length2D(edge2))] + self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, + self._edge] + elif(self.dim == 3): + xyz = self.gridN + A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, + self.nNy, + self.nNz])) + edge1 = xyz[D, :] - xyz[A, :] + A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, + self.nCy, + self.nNz])) + edge2 = xyz[B, :] - xyz[A, :] + A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, + self.nNy, + self.nCz])) + edge3 = xyz[E, :] - xyz[A, :] + self._edge = np.r_[Utils.mkvc(length3D(edge1)), + Utils.mkvc(length3D(edge2)), + Utils.mkvc(length3D(edge3))] + self._tangents = (np.r_[edge1, edge2, edge3] / + np.c_[self._edge, self._edge, self._edge]) return self._edge - return locals() - _edge = None - edge = property(**edge()) + return self._edge - def tangents(): - doc = "Edge tangents." - - def fget(self): - if(self._tangents is None): - self.edge # calling .edge will create the tangents - return self._tangents - return locals() - _tangents = None - tangents = property(**tangents()) + @property + def tangents(self): + """ + Edge tangents + """ + if getattr(self, '_tangents', None) is None: + self.edge # calling .edge will create the tangents + return self._tangents diff --git a/tests/mesh/test_CurvilinearMesh.py b/tests/mesh/test_CurvilinearMesh.py index 42e3d877..a2d15708 100644 --- a/tests/mesh/test_CurvilinearMesh.py +++ b/tests/mesh/test_CurvilinearMesh.py @@ -10,7 +10,9 @@ class BasicCurvTests(unittest.TestCase): a = np.array([1, 1, 1]) b = np.array([1, 2]) c = np.array([1, 4]) - gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h] + + def gridIt(h): return [np.cumsum(np.r_[0, x]) for x in h] + X, Y = ndgrid(gridIt([a, b]), vector=False) self.TM2 = TensorMesh([a, b]) self.Curv2 = CurvilinearMesh([X, Y]) @@ -19,7 +21,10 @@ class BasicCurvTests(unittest.TestCase): self.Curv3 = CurvilinearMesh([X, Y, Z]) 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]) self.assertTrue(np.all(self.Curv3.area == test_area)) def test_vol_3D(self): @@ -33,54 +38,85 @@ class BasicCurvTests(unittest.TestCase): 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.Curv3.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.Curv2.edge == test_edge) self.assertTrue(t1) def test_tangents(self): T = self.Curv2.tangents - self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx))) - self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx))) - self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy))) - self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == + np.ones(self.Curv2.nEx))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == + np.zeros(self.Curv2.nEx))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == + np.zeros(self.Curv2.nEy))) + self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == + np.ones(self.Curv2.nEy))) T = self.Curv3.tangents - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx))) - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx))) - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == + np.ones(self.Curv3.nEx))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == + np.zeros(self.Curv3.nEx))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == + np.zeros(self.Curv3.nEx))) - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy))) - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy))) - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == + np.zeros(self.Curv3.nEy))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == + np.ones(self.Curv3.nEy))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == + np.zeros(self.Curv3.nEy))) - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz))) - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz))) - self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == + np.zeros(self.Curv3.nEz))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == + np.zeros(self.Curv3.nEz))) + self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == + np.ones(self.Curv3.nEz))) def test_normals(self): N = self.Curv2.normals - self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx))) - self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx))) - self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy))) - self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == + np.ones(self.Curv2.nFx))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == + np.zeros(self.Curv2.nFx))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == + np.zeros(self.Curv2.nFy))) + self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == + np.ones(self.Curv2.nFy))) N = self.Curv3.normals - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx))) - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx))) - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == + np.ones(self.Curv3.nFx))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == + np.zeros(self.Curv3.nFx))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == + np.zeros(self.Curv3.nFx))) - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy))) - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy))) - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == + np.zeros(self.Curv3.nFy))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == + np.ones(self.Curv3.nFy))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == + np.zeros(self.Curv3.nFy))) - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz))) - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz))) - self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == + np.zeros(self.Curv3.nFz))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == + np.zeros(self.Curv3.nFz))) + self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == + np.ones(self.Curv3.nFz))) def test_grid(self): self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC))