From 2bdc3958b5a54bb3d3dd525b5d937cf6a8c1d120 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 16 Feb 2014 14:42:43 -0800 Subject: [PATCH] Split out base tensor mesh stuff into different mixing class. --- SimPEG/Mesh/CylMesh.py | 57 +++++--- SimPEG/Mesh/TensorMesh.py | 252 +++++++++++++++++---------------- SimPEG/Tests/test_operators.py | 1 + 3 files changed, 169 insertions(+), 141 deletions(-) diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index 34f19c1e..eb557735 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -2,35 +2,21 @@ import numpy as np import scipy.sparse as sp from scipy.constants import pi from SimPEG.Utils import mkvc, ndgrid, sdiag -from TensorMesh import TensorMesh +from TensorMesh import BaseTensorMesh -class CylMesh(TensorMesh): +class CylMesh(BaseTensorMesh): """ CylMesh is a mesh class for cylindrical problems """ _meshType = 'CYL' + _unitDimensions = [1, 2*np.pi, 1] + def __init__(self, h, x0=None): - assert len(h) == 3, "len(h) must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]" - - if x0 is not None: - assert type(x0) == np.ndarray, "x0 must be an ndarray" - assert x0.size == 3, "x0 must have 3 elements" - else: - x0 = np.r_[0, 0, 0] - - for i, h_i in enumerate(h): - if type(h_i) in [int, long, float]: - # This gives you something over the unit cylinder. - h_i = (2*np.pi if i==1 else 1.)*np.ones(int(h_i))/int(h_i) - assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) - assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) - h[i] = h_i[:] # make a copy. - - assert h[1].sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" - - TensorMesh.__init__(self, h, x0) + BaseTensorMesh.__init__(self, h, x0) + assert self.dim == 3, "dim of mesh must equal 3, for a cylindrically symmetric mesh use [hx, 1, hz]" + assert self.hy.sum() == 2*np.pi, "The 2nd dimension must sum to 2*pi" @property def nNx(self): @@ -185,6 +171,35 @@ class CylMesh(TensorMesh): # Operators #################################################### + @property + def faceDiv(self): + """Construct divergence operator (face-stg to cell-centres).""" + raise NotImplementedError('faceDiv not yet implemented') + @property + def faceDivx(self): + """Construct divergence operator in the x component (face-stg to cell-centres).""" + raise NotImplementedError('faceDivx not yet implemented') + + @property + def faceDivy(self): + """Construct divergence operator in the y component (face-stg to cell-centres).""" + raise NotImplementedError('faceDivy not yet implemented') + + @property + def faceDivz(self): + """Construct divergence operator in the z component (face-stg to cell-centres).""" + raise NotImplementedError('faceDivz not yet implemented') + + @property + def nodalGrad(self): + """Construct gradient operator (nodes to edges).""" + raise NotImplementedError('nodalGrad not yet implemented') + + @property + def nodalLaplacian(self): + """Construct laplacian operator (nodes to edges).""" + raise NotImplementedError('nodalLaplacian not yet implemented') + @property def edgeCurl(self): """The edgeCurl property.""" diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 693359bf..6f292d0f 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -4,38 +4,13 @@ from TensorView import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts -class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): - """ - TensorMesh is a mesh class that deals with tensor product meshes. - - Any Mesh that has a constant width along the entire axis - such that it can defined by a single width vector, called 'h'. - - :: - - hx = np.array([1,1,1]) - hy = np.array([1,2]) - hz = np.array([1,1,1,1]) - - mesh = Mesh.TensorMesh([hx, hy, hz]) - - Example of a padded tensor mesh: - - .. plot:: - - from SimPEG import Mesh, Utils - M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) - M.plotGrid() - - For a quick tensor mesh on a (10x12x15) unit cube:: - - mesh = Mesh.TensorMesh([10, 12, 15]) - - """ +class BaseTensorMesh(BaseRectangularMesh): __metaclass__ = Utils.SimPEGMetaClass - _meshType = 'TENSOR' + _meshType = 'BASETENSOR' + + _unitDimensions = [1, 1, 1] def __init__(self, h_in, x0=None): assert type(h_in) is list, 'h_in must be a list' @@ -43,7 +18,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): for i, h_i in enumerate(h_in): if type(h_i) in [int, long, float]: # This gives you something over the unit cube. - h_i = np.ones(int(h_i))/int(h_i) + h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i) assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) h[i] = h_i[:] # make a copy. @@ -54,57 +29,6 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): # Ensure h contains 1D vectors self._h = [Utils.mkvc(x.astype(float)) for x in h] - def __str__(self): - outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) - def printH(hx, outStr=''): - i = -1 - while True: - i = i + 1 - if i > hx.size: - break - elif i == hx.size: - break - h = hx[i] - n = 1 - for j in range(i+1, hx.size): - if hx[j] == h: - n = n + 1 - i = i + 1 - else: - break - - if n == 1: - outStr = outStr + ' {0:.2f},'.format(h) - else: - outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h) - - return outStr[:-1] - - if self.dim == 1: - outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) - outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) - outStr = outStr + printH(self.hx, outStr='\n hx:') - pass - elif self.dim == 2: - outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) - outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1]) - outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) - outStr = outStr + '\n nCy: {0:d}'.format(self.nCy) - outStr = outStr + printH(self.hx, outStr='\n hx:') - outStr = outStr + printH(self.hy, outStr='\n hy:') - elif self.dim == 3: - outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) - outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1]) - outStr = outStr + '\n z0: {0:.2f}'.format(self.x0[2]) - outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) - outStr = outStr + '\n nCy: {0:d}'.format(self.nCy) - outStr = outStr + '\n nCz: {0:d}'.format(self.nCz) - outStr = outStr + printH(self.hx, outStr='\n hx:') - outStr = outStr + printH(self.hy, outStr='\n hy:') - outStr = outStr + printH(self.hz, outStr='\n hz:') - - return outStr - @property def h(self): """h is a list containing the cell widths of the tensor mesh in each dimension.""" @@ -211,6 +135,133 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): self._gridEz = Utils.ndgrid(self.getTensor('Ez')) return self._gridEz + def getTensor(self, locType): + """ Returns a tensor list. + + :param str locType: What tensor (see below) + :rtype: list + :return: list of the tensors that make up the mesh. + + locType can be:: + + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'N' -> scalar field defined on nodes + 'CC' -> scalar field defined on cell centers + """ + + if locType is 'Fx': + ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + elif locType is 'Fy': + ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + elif locType is 'Fz': + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + elif locType is 'Ex': + ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + elif locType is 'Ey': + ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + elif locType is 'Ez': + ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + elif locType is 'CC': + ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] + elif locType is 'N': + ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + + return [t for t in ten if t is not None] + + +class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): + """ + TensorMesh is a mesh class that deals with tensor product meshes. + + Any Mesh that has a constant width along the entire axis + such that it can defined by a single width vector, called 'h'. + + :: + + hx = np.array([1,1,1]) + hy = np.array([1,2]) + hz = np.array([1,1,1,1]) + + mesh = Mesh.TensorMesh([hx, hy, hz]) + + Example of a padded tensor mesh: + + .. plot:: + + from SimPEG import Mesh, Utils + M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0)))) + M.plotGrid() + + For a quick tensor mesh on a (10x12x15) unit cube:: + + mesh = Mesh.TensorMesh([10, 12, 15]) + + """ + + __metaclass__ = Utils.SimPEGMetaClass + + _meshType = 'TENSOR' + + def __init__(self, h_in, x0=None): + BaseTensorMesh.__init__(self, h_in, x0) + + def __str__(self): + outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) + def printH(hx, outStr=''): + i = -1 + while True: + i = i + 1 + if i > hx.size: + break + elif i == hx.size: + break + h = hx[i] + n = 1 + for j in range(i+1, hx.size): + if hx[j] == h: + n = n + 1 + i = i + 1 + else: + break + + if n == 1: + outStr = outStr + ' {0:.2f},'.format(h) + else: + outStr = outStr + ' {0:d}*{1:.2f},'.format(n,h) + + return outStr[:-1] + + if self.dim == 1: + outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) + outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) + outStr = outStr + printH(self.hx, outStr='\n hx:') + pass + elif self.dim == 2: + outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) + outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1]) + outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) + outStr = outStr + '\n nCy: {0:d}'.format(self.nCy) + outStr = outStr + printH(self.hx, outStr='\n hx:') + outStr = outStr + printH(self.hy, outStr='\n hy:') + elif self.dim == 3: + outStr = outStr + '\n x0: {0:.2f}'.format(self.x0[0]) + outStr = outStr + '\n y0: {0:.2f}'.format(self.x0[1]) + outStr = outStr + '\n z0: {0:.2f}'.format(self.x0[2]) + outStr = outStr + '\n nCx: {0:d}'.format(self.nCx) + outStr = outStr + '\n nCy: {0:d}'.format(self.nCy) + outStr = outStr + '\n nCz: {0:d}'.format(self.nCz) + outStr = outStr + printH(self.hx, outStr='\n hx:') + outStr = outStr + printH(self.hy, outStr='\n hy:') + outStr = outStr + printH(self.hz, outStr='\n hz:') + + return outStr + + # --------------- Geometries --------------------- @property def vol(self): @@ -274,45 +325,6 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): # --------------- Methods --------------------- - def getTensor(self, locType): - """ Returns a tensor list. - - :param str locType: What tensor (see below) - :rtype: list - :return: list of the tensors that make up the mesh. - - locType can be:: - - 'Ex' -> x-component of field defined on edges - 'Ey' -> y-component of field defined on edges - 'Ez' -> z-component of field defined on edges - 'Fx' -> x-component of field defined on faces - 'Fy' -> y-component of field defined on faces - 'Fz' -> z-component of field defined on faces - 'N' -> scalar field defined on nodes - 'CC' -> scalar field defined on cell centers - """ - - if locType is 'Fx': - ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] - elif locType is 'Fy': - ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] - elif locType is 'Fz': - ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] - elif locType is 'Ex': - ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] - elif locType is 'Ey': - ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] - elif locType is 'Ez': - ten = [self.vectorNx , self.vectorNy , self.vectorCCz] - elif locType is 'CC': - ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] - elif locType is 'N': - ten = [self.vectorNx , self.vectorNy , self.vectorNz ] - - return [t for t in ten if t is not None] - - def isInside(self, pts): """ Determines if a set of points are inside a mesh. diff --git a/SimPEG/Tests/test_operators.py b/SimPEG/Tests/test_operators.py index 67741aa7..76884b96 100644 --- a/SimPEG/Tests/test_operators.py +++ b/SimPEG/Tests/test_operators.py @@ -3,6 +3,7 @@ import unittest from TestUtils import OrderTest import matplotlib.pyplot as plt +#TODO: 'randomTensorMesh' MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])