mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 13:13:14 +08:00
use @property decorator in curve mesh (see: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d?groups=code_patterns%3A3bECxdfc%3Af0&tab=basics)
This commit is contained in:
Lindsey Heagy
+234
-201
@@ -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
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -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))
|
||||
|
||||
Reference in New Issue
Block a user