diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 9e7c1d93..43d6c31f 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -78,29 +78,208 @@ class InnerProducts(object): def __init__(self): raise Exception('InnerProducts is a base class providing inner product matrices for meshes and cannot run on its own. Inherit to your favorite Mesh class.') - def getFaceInnerProduct(self, mu=None, returnP=False): - """Wrapper function, - - :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct` - - :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D` + def getFaceInnerProduct(M, mu=None, returnP=False): """ - if self.dim == 2: - return getFaceInnerProduct2D(self, mu, returnP) - elif self.dim == 3: - return getFaceInnerProduct(self, mu, returnP) + :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (sum(nF), sum(nF)) - def getEdgeInnerProduct(self, sigma=None, returnP=False): - """Wrapper function, + Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: - :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct` + .. math:: + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1} \end{matrix}\\right] + + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3} \end{matrix}\\right] + + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\\right] + + \mathbf{M}(\\vec{\mu}) = {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P000, P100, P010, P110, P001, P101, P011, P111] + + Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + **For 2D:** + + Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows: + + .. math:: + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right] + + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right] + + \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right] + + + .. math:: + + \mathbf{M}(\\vec{\mu}) = {1\over 4} + \left(\sum_{i=1}^4 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P00, P10, P01, P11] + + Here each P (2*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. - :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D` """ - if self.dim == 2: - return getEdgeInnerProduct2D(self, sigma, returnP) - elif self.dim == 3: - return getEdgeInnerProduct(self, sigma, returnP) + if M.dim == 2: + # Square root of cell volume multiplied by 1/4 + v = np.sqrt(0.25*M.vol) + V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry + + Pxx = _getFacePxx(M) + P000 = V2*Pxx('fXm', 'fYm') + P100 = V2*Pxx('fXp', 'fYm') + P010 = V2*Pxx('fXm', 'fYp') + P110 = V2*Pxx('fXp', 'fYp') + elif M.dim == 3: + # Square root of cell volume multiplied by 1/8 + v = np.sqrt(0.125*M.vol) + V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry + + Pxxx = _getFacePxxx(M) + P000 = V3*Pxxx('fXm', 'fYm', 'fZm') + P100 = V3*Pxxx('fXp', 'fYm', 'fZm') + P010 = V3*Pxxx('fXm', 'fYp', 'fZm') + P110 = V3*Pxxx('fXp', 'fYp', 'fZm') + P001 = V3*Pxxx('fXm', 'fYm', 'fZp') + P101 = V3*Pxxx('fXp', 'fYm', 'fZp') + P011 = V3*Pxxx('fXm', 'fYp', 'fZp') + P111 = V3*Pxxx('fXp', 'fYp', 'fZp') + + Mu = _makeTensor(M, mu) + A = P000.T*Mu*P000 + P100.T*Mu*P100 + P010.T*Mu*P010 + P110.T*Mu*P110 + P = [P000, P100, P010, P110] + if M.dim == 3: + A = A + P001.T*Mu*P001 + P101.T*Mu*P101 + P011.T*Mu*P011 + P111.T*Mu*P111 + P += [P001, P101, P011, P111] + if returnP: + return A, P + else: + return A + + def getEdgeInnerProduct(M, sigma=None, returnP=False): + """ + :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param bool returnP: returns the projection matrices + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (sum(nE), sum(nE)) + + + Depending on the number of columns (either 1, 3, or 6) of sigma, the material property is interpreted as follows: + + .. math:: + \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1} \end{matrix}\\right] + + \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3} \end{matrix}\\right] + + \Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3} \end{matrix}\\right] + + What is returned: + + .. math:: + \mathbf{M}(\Sigma) = {1\over 8} + \left(\sum_{i=1}^8 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P000, P100, P010, P110, P001, P101, P011, P111] + + Here each P (3*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + **For 2D:** + + Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows: + + .. math:: + \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right] + + \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right] + + \Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right] + + + .. math:: + + \mathbf{M}(\Sigma) = {1\over 4} + \left(\sum_{i=1}^4 + \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c + \\right) + + + If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + + P = [P00, P10, P01, P11] + + Here each P (2*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: + + .. math:: + \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} + + Note that this is completed for each cell in the mesh at the same time. + + """ + # We will multiply by V on each side to keep symmetry + if M.dim == 2: + # Square root of cell volume multiplied by 1/4 + v = np.sqrt(0.25*M.vol) + V = sdiag(np.r_[v, v]) + eP = _getEdgePxx(M) + P000 = V*eP('eX0', 'eY0') + P100 = V*eP('eX0', 'eY1') + P010 = V*eP('eX1', 'eY0') + P110 = V*eP('eX1', 'eY1') + elif M.dim == 3: + # Square root of cell volume multiplied by 1/8 + v = np.sqrt(0.125*M.vol) + V = sdiag(np.r_[v, v, v]) + eP = _getEdgePxxx(M) + P000 = V*eP('eX0', 'eY0', 'eZ0') + P100 = V*eP('eX0', 'eY1', 'eZ1') + P010 = V*eP('eX1', 'eY0', 'eZ2') + P110 = V*eP('eX1', 'eY1', 'eZ3') + P001 = V*eP('eX2', 'eY2', 'eZ0') + P101 = V*eP('eX2', 'eY3', 'eZ1') + P011 = V*eP('eX3', 'eY2', 'eZ2') + P111 = V*eP('eX3', 'eY3', 'eZ3') + + Sigma = _makeTensor(M, sigma) + A = P000.T*Sigma*P000 + P100.T*Sigma*P100 + P010.T*Sigma*P010 + P110.T*Sigma*P110 + P = [P000, P100, P010, P110] + if M.dim == 3: + A = A + P001.T*Sigma*P001 + P101.T*Sigma*P101 + P011.T*Sigma*P011 + P111.T*Sigma*P111 + P += [P001, P101, P011, P111] + if returnP: + return A, P + else: + return A # ------------------------ Geometries ------------------------------ # @@ -121,434 +300,264 @@ class InnerProducts(object): # | |/ # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) +def _makeTensor(M, sigma): + if sigma is None: # default is ones + sigma = np.ones((M.nC, 1)) -def getFaceInnerProduct(mesh, mu=None, returnP=False): - """ - :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (sum(nF), sum(nF)) + if M.dim == 2: + if sigma.size == M.nC: # Isotropic! + sigma = mkvc(sigma) # ensure it is a vector. + Sigma = sdiag(np.r_[sigma, sigma]) + elif sigma.shape[1] == 2: # Diagonal tensor + Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]]) + elif sigma.shape[1] == 3: # Fully anisotropic + row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2]))) + row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1]))) + Sigma = sp.vstack((row1, row2)) + elif M.dim == 3: + if sigma.size == M.nC: # Isotropic! + sigma = mkvc(sigma) # ensure it is a vector. + Sigma = sdiag(np.r_[sigma, sigma, sigma]) + elif sigma.shape[1] == 3: # Diagonal tensor + Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) + elif sigma.shape[1] == 6: # Fully anisotropic + row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) + row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) + row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) + Sigma = sp.vstack((row1, row2, row3)) + return Sigma - Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows: +def _getFacePxx(M): + if M._meshType == 'TREE': + return M._getFacePxx - .. math:: - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{1} & 0 \\\\ 0 & 0 & \mu_{1} \end{matrix}\\right] + return _getFacePxx_Rectangular(M) - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 & 0 \\\\ 0 & \mu_{2} & 0 \\\\ 0 & 0 & \mu_{3} \end{matrix}\\right] +def _getFacePxxx(M): + if M._meshType == 'TREE': + return M._getFacePxxx - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{4} & \mu_{5} \\\\ \mu_{4} & \mu_{2} & \mu_{6} \\\\ \mu_{5} & \mu_{6} & \mu_{3} \end{matrix}\\right] + return _getFacePxxx_Rectangular(M) - \mathbf{M}(\\vec{\mu}) = {1\over 8} - \left(\sum_{i=1}^8 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c - \\right) +def _getEdgePxx(M): + if M._meshType == 'TREE': + return M._getEdgePxx - If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: + return _getEdgePxx_Rectangular(M) - P = [P000, P001, P010, P011, P100, P101, P110, P111] +def _getEdgePxxx(M): + if M._meshType == 'TREE': + return M._getEdgePxxx - Here each P (3*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: + return _getEdgePxxx_Rectangular(M) - .. math:: - \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} +def _getFacePxx_Rectangular(M): + """returns a function for creating projection matrices - Note that this is completed for each cell in the mesh at the same time. + Mats takes you from faces a subset of all faces on only the + faces that you ask for. - """ + These are centered around a single nodes. - if mu is None: # default is ones - mu = np.ones((mesh.nC, 1)) + For example, if this was your entire mesh: - m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) - nc = mesh.nC + f3(Yp) + 2_______________3 + | | + | | + | | + f0(Xm) | x | f1(Xp) + | | + | | + |_______________| + 0 1 + f2(Ym) - i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) + Pxx('m','m') = | 1, 0, 0, 0 | + | 0, 0, 1, 0 | - iijjkk = ndgrid(i, j, k) - ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] + Pxx('p','m') = | 0, 1, 0, 0 | + | 0, 0, 1, 0 | - if mesh._meshType == 'LOM': - fN1 = mesh.r(mesh.normals, 'F', 'Fx', 'M') - fN2 = mesh.r(mesh.normals, 'F', 'Fy', 'M') - fN3 = mesh.r(mesh.normals, 'F', 'Fz', 'M') - - def Pxxx(pos): - ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nFv[0] - ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nFv[0] + mesh.nFv[1] - - IND = np.r_[ind1, ind2, ind3].flatten() - - PXXX = sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nF))).tocsr() - - if mesh._meshType == 'LOM': - I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(fN1[1], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(fN1[2], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), - getSubArray(fN2[0], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(fN2[1], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(fN2[2], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), - getSubArray(fN3[0], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(fN3[1], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(fN3[2], [i + pos[2][0], j + pos[2][1], k + pos[2][2]])) - PXXX = I3x3 * PXXX - - return PXXX - - # no | node | f1 | f2 | f3 - # 000 | i ,j ,k | i , j, k | i, j , k | i, j, k - # 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k - # 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k - # 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k - # 001 | i ,j ,k+1 | i , j, k | i, j , k | i, j, k+1 - # 101 | i+1,j ,k+1 | i+1, j, k | i, j , k | i, j, k+1 - # 011 | i ,j+1,k+1 | i , j, k | i, j+1, k | i, j, k+1 - # 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1 - - # Square root of cell volume multiplied by 1/8 - v = np.sqrt(0.125*mesh.vol) - V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry - - P000 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) - P100 = V3*Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 0]]) - P010 = V3*Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 0]]) - P110 = V3*Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 0]]) - P001 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 1]]) - P101 = V3*Pxxx([[1, 0, 0], [0, 0, 0], [0, 0, 1]]) - P011 = V3*Pxxx([[0, 0, 0], [0, 1, 0], [0, 0, 1]]) - P111 = V3*Pxxx([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) - - if mu.size == mesh.nC: # Isotropic! - mu = mkvc(mu) # ensure it is a vector. - Mu = sdiag(np.r_[mu, mu, mu]) - elif mu.shape[1] == 3: # Diagonal tensor - Mu = sdiag(np.r_[mu[:, 0], mu[:, 1], mu[:, 2]]) - elif mu.shape[1] == 6: # Fully anisotropic - row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 3]), sdiag(mu[:, 4]))) - row2 = sp.hstack((sdiag(mu[:, 3]), sdiag(mu[:, 1]), sdiag(mu[:, 5]))) - row3 = sp.hstack((sdiag(mu[:, 4]), sdiag(mu[:, 5]), sdiag(mu[:, 2]))) - Mu = sp.vstack((row1, row2, row3)) - - A = P000.T*Mu*P000 + P001.T*Mu*P001 + P010.T*Mu*P010 + P011.T*Mu*P011 + P100.T*Mu*P100 + P101.T*Mu*P101 + P110.T*Mu*P110 + P111.T*Mu*P111 - P = [P000, P001, P010, P011, P100, P101, P110, P111] - if returnP: - return A, P - else: - return A - - -def getFaceInnerProduct2D(mesh, mu=None, returnP=False): - """ - :param numpy.array mu: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3)) - :param bool returnP: returns the projection matrices - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (sum(nF), sum(nF)) - - Depending on the number of columns (either 1, 2, or 3) of mu, the material property is interpreted as follows: - - .. math:: - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{1} \end{matrix}\\right] - - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & 0 \\\\ 0 & \mu_{2} \end{matrix}\\right] - - \\vec{\mu} = \left[\\begin{matrix} \mu_{1} & \mu_{3} \\\\ \mu_{3} & \mu_{2} \end{matrix}\\right] - - - .. math:: - - \mathbf{M}(\\vec{\mu}) = {1\over 4} - \left(\sum_{i=1}^4 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \\vec{\mu} \sqrt{v_{\\text{cell}}} \mathbf{J}_c - \\right) - - - If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: - - P = [P00, P10, P01, P11] - - Here each P (2*nC, sum(nF)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: - - .. math:: - \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} - - Note that this is completed for each cell in the mesh at the same time. - - """ - - if mu is None: # default is ones - mu = np.ones((mesh.nC, 1)) - - m = np.array([mesh.nCx, mesh.nCy]) - nc = mesh.nC - - i, j = np.int64(range(m[0])), np.int64(range(m[1])) + """ + i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy)) iijj = ndgrid(i, j) ii, jj = iijj[:, 0], iijj[:, 1] - if mesh._meshType == 'LOM': - fN1 = mesh.r(mesh.normals, 'F', 'Fx', 'M') - fN2 = mesh.r(mesh.normals, 'F', 'Fy', 'M') + if M._meshType == 'LOM': + fN1 = M.r(M.normals, 'F', 'Fx', 'M') + fN2 = M.r(M.normals, 'F', 'Fy', 'M') - def Pxx(pos): - ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1]]) - ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nFv[0] + def Pxx(xFace, yFace): + """ + xFace is 'fXp' or 'fXm' + yFace is 'fYp' or 'fYm' + """ + # no | node | f1 | f2 + # 00 | i ,j | i , j | i, j + # 10 | i+1,j | i+1, j | i, j + # 01 | i ,j+1 | i , j | i, j+1 + # 11 | i+1,j+1 | i+1, j | i, j+1 + + posFx = 0 if xFace == 'fXm' else 1 + posFy = 0 if yFace == 'fYm' else 1 + + ind1 = sub2ind(M.nFx, np.c_[ii + posFx, jj]) + ind2 = sub2ind(M.nFy, np.c_[ii, jj + posFy]) + M.nFv[0] IND = np.r_[ind1, ind2].flatten() - PXX = sp.coo_matrix((np.ones(2*nc), (range(2*nc), IND)), shape=(2*nc, np.sum(mesh.nF))).tocsr() + PXX = sp.csr_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nF))) - if mesh._meshType == 'LOM': - I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + pos[0][0], j + pos[0][1]]), getSubArray(fN1[1], [i + pos[0][0], j + pos[0][1]]), - getSubArray(fN2[0], [i + pos[1][0], j + pos[1][1]]), getSubArray(fN2[1], [i + pos[1][0], j + pos[1][1]])) + if M._meshType == 'LOM': + I2x2 = inv2X2BlockDiagonal(getSubArray(fN1[0], [i + posFx, j]), getSubArray(fN1[1], [i + posFx, j]), + getSubArray(fN2[0], [i, j + posFy]), getSubArray(fN2[1], [i, j + posFy])) PXX = I2x2 * PXX return PXX - # no | node | f1 | f2 - # 00 | i ,j | i , j | i, j - # 10 | i+1,j | i+1, j | i, j - # 01 | i ,j+1 | i , j | i, j+1 - # 11 | i+1,j+1 | i+1, j | i, j+1 + return Pxx - # Square root of cell volume multiplied by 1/4 - v = np.sqrt(0.25*mesh.vol) - V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry +def _getFacePxxx_Rectangular(M): + """returns a function for creating projection matrices - P00 = V2*Pxx([[0, 0], [0, 0]]) - P10 = V2*Pxx([[1, 0], [0, 0]]) - P01 = V2*Pxx([[0, 0], [0, 1]]) - P11 = V2*Pxx([[1, 0], [0, 1]]) + Mats takes you from faces a subset of all faces on only the + faces that you ask for. - if mu.size == mesh.nC: # Isotropic! - mu = mkvc(mu) # ensure it is a vector. - Mu = sdiag(np.r_[mu, mu]) - elif mu.shape[1] == 2: # Diagonal tensor - Mu = sdiag(np.r_[mu[:, 0], mu[:, 1]]) - elif mu.shape[1] == 3: # Fully anisotropic - row1 = sp.hstack((sdiag(mu[:, 0]), sdiag(mu[:, 2]))) - row2 = sp.hstack((sdiag(mu[:, 2]), sdiag(mu[:, 1]))) - Mu = sp.vstack((row1, row2)) - - A = P00.T*Mu*P00 + P10.T*Mu*P10 + P01.T*Mu*P01 + P11.T*Mu*P11 - P = [P00, P10, P01, P11] - if returnP: - return A, P - else: - return A - - -def getEdgeInnerProduct(mesh, sigma=None, returnP=False): - """ - :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param bool returnP: returns the projection matrices - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (sum(nE), sum(nE)) - - - Depending on the number of columns (either 1, 3, or 6) of sigma, the material property is interpreted as follows: - - .. math:: - \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{1} & 0 \\\\ 0 & 0 & \sigma_{1} \end{matrix}\\right] - - \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 & 0 \\\\ 0 & \sigma_{2} & 0 \\\\ 0 & 0 & \sigma_{3} \end{matrix}\\right] - - \Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{4} & \sigma_{5} \\\\ \sigma_{4} & \sigma_{2} & \sigma_{6} \\\\ \sigma_{5} & \sigma_{6} & \sigma_{3} \end{matrix}\\right] - - What is returned: - - .. math:: - \mathbf{M}(\Sigma) = {1\over 8} - \left(\sum_{i=1}^8 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c - \\right) - - If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: - - P = [P000, P001, P010, P011, P100, P101, P110, P111] - - Here each P (3*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: - - .. math:: - \mathbf{P}_{(i)} = \sqrt{ {1\over 8} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} - - Note that this is completed for each cell in the mesh at the same time. + These are centered around a single nodes. """ - if sigma is None: # default is ones - sigma = np.ones((mesh.nC, 1)) - - m = np.array([mesh.nCx, mesh.nCy, mesh.nCz]) - nc = mesh.nC - - i, j, k = np.int64(range(m[0])), np.int64(range(m[1])), np.int64(range(m[2])) + i, j, k = np.int64(range(M.nCx)), np.int64(range(M.nCy)), np.int64(range(M.nCz)) iijjkk = ndgrid(i, j, k) ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - if mesh._meshType == 'LOM': - eT1 = mesh.r(mesh.tangents, 'E', 'Ex', 'M') - eT2 = mesh.r(mesh.tangents, 'E', 'Ey', 'M') - eT3 = mesh.r(mesh.tangents, 'E', 'Ez', 'M') + if M._meshType == 'LOM': + fN1 = M.r(M.normals, 'F', 'Fx', 'M') + fN2 = M.r(M.normals, 'F', 'Fy', 'M') + fN3 = M.r(M.normals, 'F', 'Fz', 'M') - def Pxxx(pos): - ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]]) - ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nEv[0] - ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nEv[0] + mesh.nEv[1] + def Pxxx(xFace, yFace, zFace): + """ + xFace is 'fXp' or 'fXm' + yFace is 'fYp' or 'fYm' + zFace is 'fZp' or 'fZm' + """ + + # no | node | f1 | f2 | f3 + # 000 | i ,j ,k | i , j, k | i, j , k | i, j, k + # 100 | i+1,j ,k | i+1, j, k | i, j , k | i, j, k + # 010 | i ,j+1,k | i , j, k | i, j+1, k | i, j, k + # 110 | i+1,j+1,k | i+1, j, k | i, j+1, k | i, j, k + # 001 | i ,j ,k+1 | i , j, k | i, j , k | i, j, k+1 + # 101 | i+1,j ,k+1 | i+1, j, k | i, j , k | i, j, k+1 + # 011 | i ,j+1,k+1 | i , j, k | i, j+1, k | i, j, k+1 + # 111 | i+1,j+1,k+1 | i+1, j, k | i, j+1, k | i, j, k+1 + + posX = 0 if xFace == 'fXm' else 1 + posY = 0 if yFace == 'fYm' else 1 + posZ = 0 if zFace == 'fZm' else 1 + + ind1 = sub2ind(M.nFx, np.c_[ii + posX, jj, kk]) + ind2 = sub2ind(M.nFy, np.c_[ii, jj + posY, kk]) + M.nFv[0] + ind3 = sub2ind(M.nFz, np.c_[ii, jj, kk + posZ]) + M.nFv[0] + M.nFv[1] IND = np.r_[ind1, ind2, ind3].flatten() - PXXX = sp.coo_matrix((np.ones(3*nc), (range(3*nc), IND)), shape=(3*nc, np.sum(mesh.nE))).tocsr() + PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nF))).tocsr() - if mesh._meshType == 'LOM': - I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(eT1[1], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), getSubArray(eT1[2], [i + pos[0][0], j + pos[0][1], k + pos[0][2]]), - getSubArray(eT2[0], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(eT2[1], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), getSubArray(eT2[2], [i + pos[1][0], j + pos[1][1], k + pos[1][2]]), - getSubArray(eT3[0], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(eT3[1], [i + pos[2][0], j + pos[2][1], k + pos[2][2]]), getSubArray(eT3[2], [i + pos[2][0], j + pos[2][1], k + pos[2][2]])) + if M._meshType == 'LOM': + I3x3 = inv3X3BlockDiagonal(getSubArray(fN1[0], [i + posX, j, k]), getSubArray(fN1[1], [i + posX, j, k]), getSubArray(fN1[2], [i + posX, j, k]), + getSubArray(fN2[0], [i, j + posY, k]), getSubArray(fN2[1], [i, j + posY, k]), getSubArray(fN2[2], [i, j + posY, k]), + getSubArray(fN3[0], [i, j, k + posZ]), getSubArray(fN3[1], [i, j, k + posZ]), getSubArray(fN3[2], [i, j, k + posZ])) PXXX = I3x3 * PXXX return PXXX + return Pxxx - # no | node | e1 | e2 | e3 - # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k - # 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k - # 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k - # 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k - # 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k - # 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k - # 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k - # 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k - - # Square root of cell volume multiplied by 1/8 - v = np.sqrt(0.125*mesh.vol) - V3 = sdiag(np.r_[v, v, v]) # We will multiply on each side to keep symmetry - - P000 = V3*Pxxx([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) - P100 = V3*Pxxx([[0, 0, 0], [1, 0, 0], [1, 0, 0]]) - P010 = V3*Pxxx([[0, 1, 0], [0, 0, 0], [0, 1, 0]]) - P110 = V3*Pxxx([[0, 1, 0], [1, 0, 0], [1, 1, 0]]) - P001 = V3*Pxxx([[0, 0, 1], [0, 0, 1], [0, 0, 0]]) - P101 = V3*Pxxx([[0, 0, 1], [1, 0, 1], [1, 0, 0]]) - P011 = V3*Pxxx([[0, 1, 1], [0, 0, 1], [0, 1, 0]]) - P111 = V3*Pxxx([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) - - if sigma.size == mesh.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(np.r_[sigma, sigma, sigma]) - elif sigma.shape[1] == 3: # Diagonal tensor - Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) - elif sigma.shape[1] == 6: # Fully anisotropic - row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) - row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) - row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) - Sigma = sp.vstack((row1, row2, row3)) - - A = P000.T*Sigma*P000 + P001.T*Sigma*P001 + P010.T*Sigma*P010 + P011.T*Sigma*P011 + P100.T*Sigma*P100 + P101.T*Sigma*P101 + P110.T*Sigma*P110 + P111.T*Sigma*P111 - P = [P000, P001, P010, P011, P100, P101, P110, P111] - if returnP: - return A, P - else: - return A - - -def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False): - """ - :param numpy.array sigma: material property (tensor properties are possible) at each cell center (nC, (1, 2, or 3)) - :param bool returnP: returns the projection matrices - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (sum(nE), sum(nE)) - - Depending on the number of columns (either 1, 2, or 3) of sigma, the material property is interpreted as follows: - - .. math:: - \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{1} \end{matrix}\\right] - - \Sigma = \left[\\begin{matrix} \sigma_{1} & 0 \\\\ 0 & \sigma_{2} \end{matrix}\\right] - - \Sigma = \left[\\begin{matrix} \sigma_{1} & \sigma_{3} \\\\ \sigma_{3} & \sigma_{2} \end{matrix}\\right] - - - .. math:: - - \mathbf{M}(\Sigma) = {1\over 4} - \left(\sum_{i=1}^4 - \mathbf{J}_c^{-\\top} \sqrt{v_{\\text{cell}}} \Sigma \sqrt{v_{\\text{cell}}} \mathbf{J}_c - \\right) - - - If requested (returnP=True) the projection matricies are returned as well (ordered by nodes):: - - P = [P00, P10, P01, P11] - - Here each P (2*nC, sum(nE)) is a combination of the projection, volume, and any normalization to Cartesian coordinates: - - .. math:: - \mathbf{P}_{(i)} = \sqrt{ {1\over 4} v_{\\text{cell}}} \overbrace{\mathbf{N}_{(i)}^{-1}}^{\\text{LOM only}} \mathbf{Q}_{(i)} - - Note that this is completed for each cell in the mesh at the same time. - - """ - - if sigma is None: # default is ones - sigma = np.ones((mesh.nC, 1)) - - m = np.array([mesh.nCx, mesh.nCy]) - nc = mesh.nC - - i, j = np.int64(range(m[0])), np.int64(range(m[1])) +def _getEdgePxx_Rectangular(M): + i, j = np.int64(range(M.nCx)), np.int64(range(M.nCy)) iijj = ndgrid(i, j) ii, jj = iijj[:, 0], iijj[:, 1] - if mesh._meshType == 'LOM': - eT1 = mesh.r(mesh.tangents, 'E', 'Ex', 'M') - eT2 = mesh.r(mesh.tangents, 'E', 'Ey', 'M') + if M._meshType == 'LOM': + eT1 = M.r(M.tangents, 'E', 'Ex', 'M') + eT2 = M.r(M.tangents, 'E', 'Ey', 'M') - def Pxx(pos): - ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1]]) - ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nEv[0] + def Pxx(xEdge, yEdge): + # no | node | e1 | e2 + # 00 | i ,j | i ,j | i ,j + # 10 | i+1,j | i ,j | i+1,j + # 01 | i ,j+1 | i ,j+1 | i ,j + # 11 | i+1,j+1 | i ,j+1 | i+1,j + posX = 0 if xEdge == 'eX0' else 1 + posY = 0 if yEdge == 'eY0' else 1 + + ind1 = sub2ind(M.nEx, np.c_[ii, jj + posX]) + ind2 = sub2ind(M.nEy, np.c_[ii + posY, jj]) + M.nEv[0] IND = np.r_[ind1, ind2].flatten() - PXX = sp.coo_matrix((np.ones(2*nc), (range(2*nc), IND)), shape=(2*nc, np.sum(mesh.nE))).tocsr() + PXX = sp.coo_matrix((np.ones(2*M.nC), (range(2*M.nC), IND)), shape=(2*M.nC, np.sum(M.nE))).tocsr() - if mesh._meshType == 'LOM': - I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i + pos[0][0], j + pos[0][1]]), getSubArray(eT1[1], [i + pos[0][0], j + pos[0][1]]), - getSubArray(eT2[0], [i + pos[1][0], j + pos[1][1]]), getSubArray(eT2[1], [i + pos[1][0], j + pos[1][1]])) + if M._meshType == 'LOM': + I2x2 = inv2X2BlockDiagonal(getSubArray(eT1[0], [i, j + posX]), getSubArray(eT1[1], [i, j + posX]), + getSubArray(eT2[0], [i + posY, j]), getSubArray(eT2[1], [i + posY, j])) PXX = I2x2 * PXX return PXX + return Pxx - # no | node | e1 | e2 - # 00 | i ,j | i ,j | i ,j - # 10 | i+1,j | i ,j | i+1,j - # 01 | i ,j+1 | i ,j+1 | i ,j - # 11 | i+1,j+1 | i ,j+1 | i+1,j +def _getEdgePxxx_Rectangular(M): + i, j, k = np.int64(range(M.nCx)), np.int64(range(M.nCy)), np.int64(range(M.nCz)) - # Square root of cell volume multiplied by 1/4 - v = np.sqrt(0.25*mesh.vol) - V2 = sdiag(np.r_[v, v]) # We will multiply on each side to keep symmetry + iijjkk = ndgrid(i, j, k) + ii, jj, kk = iijjkk[:, 0], iijjkk[:, 1], iijjkk[:, 2] - P00 = V2*Pxx([[0, 0], [0, 0]]) - P10 = V2*Pxx([[0, 0], [1, 0]]) - P01 = V2*Pxx([[0, 1], [0, 0]]) - P11 = V2*Pxx([[0, 1], [1, 0]]) + if M._meshType == 'LOM': + eT1 = M.r(M.tangents, 'E', 'Ex', 'M') + eT2 = M.r(M.tangents, 'E', 'Ey', 'M') + eT3 = M.r(M.tangents, 'E', 'Ez', 'M') - if sigma.size == mesh.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(np.r_[sigma, sigma]) - elif sigma.shape[1] == 2: # Diagonal tensor - Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]]) - elif sigma.shape[1] == 3: # Fully anisotropic - row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2]))) - row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1]))) - Sigma = sp.vstack((row1, row2)) + def Pxxx(xEdge, yEdge, zEdge): - A = P00.T*Sigma*P00 + P10.T*Sigma*P10 + P01.T*Sigma*P01 + P11.T*Sigma*P11 - P = [P00, P10, P01, P11] - if returnP: - return A, P - else: - return A + # no | node | e1 | e2 | e3 + # 000 | i ,j ,k | i ,j ,k | i ,j ,k | i ,j ,k + # 100 | i+1,j ,k | i ,j ,k | i+1,j ,k | i+1,j ,k + # 010 | i ,j+1,k | i ,j+1,k | i ,j ,k | i ,j+1,k + # 110 | i+1,j+1,k | i ,j+1,k | i+1,j ,k | i+1,j+1,k + # 001 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k+1 | i ,j ,k + # 101 | i+1,j ,k+1 | i ,j ,k+1 | i+1,j ,k+1 | i+1,j ,k + # 011 | i ,j+1,k+1 | i ,j+1,k+1 | i ,j ,k+1 | i ,j+1,k + # 111 | i+1,j+1,k+1 | i ,j+1,k+1 | i+1,j ,k+1 | i+1,j+1,k + posX = [0,0] if xEdge == 'eX0' else [1, 0] if xEdge == 'eX1' else [0,1] if xEdge == 'eX2' else [1,1] + posY = [0,0] if yEdge == 'eY0' else [1, 0] if yEdge == 'eY1' else [0,1] if yEdge == 'eY2' else [1,1] + posZ = [0,0] if zEdge == 'eZ0' else [1, 0] if zEdge == 'eZ1' else [0,1] if zEdge == 'eZ2' else [1,1] + + ind1 = sub2ind(M.nEx, np.c_[ii, jj + posX[0], kk + posX[1]]) + ind2 = sub2ind(M.nEy, np.c_[ii + posY[0], jj, kk + posY[1]]) + M.nEv[0] + ind3 = sub2ind(M.nEz, np.c_[ii + posZ[0], jj + posZ[1], kk]) + M.nEv[0] + M.nEv[1] + + IND = np.r_[ind1, ind2, ind3].flatten() + + PXXX = sp.coo_matrix((np.ones(3*M.nC), (range(3*M.nC), IND)), shape=(3*M.nC, np.sum(M.nE))).tocsr() + + if M._meshType == 'LOM': + I3x3 = inv3X3BlockDiagonal(getSubArray(eT1[0], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[1], [i, j + posX[0], k + posX[1]]), getSubArray(eT1[2], [i, j + posX[0], k + posX[1]]), + getSubArray(eT2[0], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[1], [i + posY[0], j, k + posY[1]]), getSubArray(eT2[2], [i + posY[0], j, k + posY[1]]), + getSubArray(eT3[0], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[1], [i + posZ[0], j + posZ[1], k]), getSubArray(eT3[2], [i + posZ[0], j + posZ[1], k])) + PXXX = I3x3 * PXXX + + return PXXX + return Pxxx if __name__ == '__main__': from TensorMesh import TensorMesh h = [np.array([1, 2, 3, 4]), np.array([1, 2, 1, 4, 2]), np.array([1, 1, 4, 1])] - mesh = TensorMesh(h) - mu = np.ones((mesh.nC, 6)) - A, P = mesh.getFaceInnerProduct(mu, returnP=True) - B, P = mesh.getEdgeInnerProduct(mu, returnP=True) + M = TensorMesh(h) + mu = np.ones((M.nC, 6)) + A, P = M.getFaceInnerProduct(mu, returnP=True) + B, P = M.getEdgeInnerProduct(mu, returnP=True) diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py new file mode 100644 index 00000000..3d2c2794 --- /dev/null +++ b/SimPEG/Mesh/TreeMesh.py @@ -0,0 +1,1116 @@ +from SimPEG import np, sp, Utils, Solver +from BaseMesh import BaseMesh +from InnerProducts import InnerProducts +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +import matplotlib.colors as colors +import matplotlib.cm as cmx + + + +def SortByX0(): + eps = 1e-7 + def mycmp(c1,c2): + if c1.x0.size == 2: + if np.abs(c1.x0[1] - c2.x0[1]) < eps: + return c1.x0[0] - c2.x0[0] + return c1.x0[1] - c2.x0[1] + elif c1.x0.size == 3: + if np.abs(c1.x0[2] - c2.x0[2]) < eps: + if np.abs(c1.x0[1] - c2.x0[1]) < eps: + return c1.x0[0] - c2.x0[0] + return c1.x0[1] - c2.x0[1] + return c1.x0[2] - c2.x0[2] + + class K(object): + def __init__(self, obj, *args): + self.obj = obj + def __lt__(self, other): + return mycmp(self.obj, other.obj) < 0 + def __gt__(self, other): + return mycmp(self.obj, other.obj) > 0 + def __eq__(self, other): + return mycmp(self.obj, other.obj) == 0 + def __le__(self, other): + return mycmp(self.obj, other.obj) <= 0 + def __ge__(self, other): + return mycmp(self.obj, other.obj) >= 0 + def __ne__(self, other): + return mycmp(self.obj, other.obj) != 0 + return K + + +class TreeNode(object): + """docstring for TreeNode""" + + __slots__ = ['x0', 'num'] + + def __init__(self, mesh, x0=[0,0]): + self.x0 = np.array(x0, dtype=float) + mesh.nodes.add(self) + + @property + def center(self): return self.x0 + +class TreeEdge(object): + """docstring for TreeEdge""" + + __slots__ = ['mesh', 'children', 'depth', 'x0', 'num', 'edgeType', 'sz', 'node0', 'node1'] + + def __init__(self, mesh, x0=[0,0], edgeType=None, sz=[1,], depth=0, + node0=None, node1=None): + self.mesh = mesh + self.depth = depth + + self.x0 = x0 + self.sz = sz + self.edgeType = edgeType + + mesh.edges.add(self) + if edgeType is 'x': mesh.edgesX.add(self) + elif edgeType is 'y': mesh.edgesY.add(self) + elif edgeType is 'z': mesh.edgesZ.add(self) + + self.node0 = node0 if isinstance(node0,TreeNode) else TreeNode(mesh, x0=self.x0) + self.node1 = node1 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=self.x0 + self.tangent*self.sz[0]) + + @property + def isleaf(self): return getattr(self, 'children', None) is None + + def refine(self): + if not self.isleaf: return + self.mesh.isNumbered = False + + self.children = np.empty(2,dtype=TreeFace) + # Create refined x0's + x0r_0 = self.x0 + x0r_1 = self.x0+0.5*self.tangent*self.sz + self.children[0] = TreeEdge(self.mesh, x0=x0r_0, edgeType=self.edgeType, sz=0.5*self.sz, depth=self.depth+1, node0=self.node0) + self.children[1] = TreeEdge(self.mesh, x0=x0r_1, edgeType=self.edgeType, sz=0.5*self.sz, depth=self.depth+1, node0=self.children[0].node1, node1=self.node1) + self.mesh.edges.remove(self) + if self.edgeType is 'x': + self.mesh.edgesX.remove(self) + elif self.edgeType is 'y': + self.mesh.edgesY.remove(self) + elif self.edgeType is 'z': + self.mesh.edgesZ.remove(self) + + @property + def tangent(self): + if self.edgeType is 'x': return np.r_[1.,0,0] + elif self.edgeType is 'y': return np.r_[0,1.,0] + elif self.edgeType is 'z': return np.r_[0,0,1.] + + def plotGrid(self, ax, text=False, lineOpts={'color':'r', 'ls': '-'}): + line = np.c_[self.node0.x0, self.node1.x0].T + ax.plot(line[:,0], line[:,1], zs=line[:,2], **lineOpts) + + @property + def center(self): + return 0.5*(self.node0.x0 + self.node1.x0) + + @property + def length(self): + return np.sqrt(((self.node1.x0 - self.node0.x0)**2).sum()) + + @property + def index(self): + if self.isleaf: return [self.num] + l = [edge.index for edge in self.children.flatten(order='F')] + # Flatten the list + # e.g. + # [[1,3],[4]] --> [1, 3, 4] + return [item for sublist in l for item in sublist] + +class TreeFace(object): + """docstring for TreeFace""" + + __slots__ = ['mesh', 'children', 'depth', 'num', 'faceType', 'sz', 'node0', 'node1', 'node2', 'node3', 'edge0', 'edge1', 'edge2', 'edge3', '_tangent0', '_tangent1'] + + def __init__(self, mesh, x0=[0,0], faceType=None, sz=[1,], depth=0, + node0=None, node1=None, + edge0=None, edge1=None, edge2=None, edge3=None): + + self.mesh = mesh + self.depth = depth + + self.faceType = faceType + self.sz = sz + + mesh.faces.add(self) + if faceType is 'x': mesh.facesX.add(self) + elif faceType is 'y': mesh.facesY.add(self) + elif faceType is 'z': mesh.facesZ.add(self) + if self.dim == 2: + # Add the nodes: + self.node0 = node0 if isinstance(node0,TreeNode) else TreeNode(mesh, x0=x0) + self.node1 = node1 if isinstance(node1,TreeNode) else TreeNode(mesh, x0=x0 + self.tangent0*self.sz[0]) + if self.dim == 3: + #TODO: Change this to edges + + # + # 2___________3 + # | e1 | + # | | + # e2| x |e3 t1 + # | | ^ + # |___________| |___> t0 + # 0 e0 1 + # + + N = {} + n0 = getattr(edge0, 'node0', None) or getattr(edge2, 'node0', None) + n1 = getattr(edge0, 'node1', None) or getattr(edge3, 'node0', None) + n2 = getattr(edge1, 'node0', None) or getattr(edge2, 'node1', None) + n3 = getattr(edge1, 'node1', None) or getattr(edge3, 'node1', None) + + eType = ['x', 'y'] if self.faceType == 'z' else ['x', 'z'] if self.faceType == 'y' else ['y', 'z'] + + e0 = edge0 if isinstance(edge0,TreeEdge) else TreeEdge(mesh, x0=x0, edgeType=eType[0], sz=np.r_[sz[0]], depth=depth, node0=n0, node1=n1) + n0, n1 = e0.node0, e0.node1 + + e1 = edge1 if isinstance(edge1,TreeEdge) else TreeEdge(mesh, x0=x0 + self.tangent1*self.sz[1], edgeType=eType[0], sz=np.r_[sz[0]], depth=depth, node0=n2, node1=n3) + n2, n3 = e1.node0, e1.node1 + + e2 = edge2 if isinstance(edge2,TreeEdge) else TreeEdge(mesh, x0=x0, edgeType=eType[1], sz=np.r_[sz[1]], depth=depth, node0=n0, node1=n2) + n0, n2 = e2.node0, e2.node1 + + e3 = edge3 if isinstance(edge3,TreeEdge) else TreeEdge(mesh, x0=x0 + self.tangent0*self.sz[0], edgeType=eType[1], sz=np.r_[sz[1]], depth=depth, node0=n1, node1=n3) + n1, n3 = e3.node0, e3.node1 + + # self.nodes = N + self.node0, self.node1, self.node2, self.node3 = n0, n1, n2, n3 + self.edge0, self.edge1, self.edge2, self.edge3 = e0, e1, e2, e3 + # self.edges = {'e0':e0, 'e1':e1, 'e2':e2, 'e3':e3} + + @property + def dim(self): return self.mesh.dim + + @property + def x0(self): return self.node0.x0 + + @property + def isleaf(self): return getattr(self, 'children', None) is None + + @property + def branchdepth(self): + if self.isleaf: + return self.depth + else: + return np.max([node.branchdepth for node in self.children.flatten('F')]) + + @property + def tangent0(self): + if getattr(self,'_tangent0',None) is None: + if self.faceType is 'x': t = np.r_[0,1.,0] + elif self.faceType is 'y': t = np.r_[1.,0,0] + elif self.faceType is 'z': t = np.r_[1.,0,0] + self._tangent0 = t[:self.dim] + return self._tangent0 + + @property + def tangent1(self): + if self.dim == 2: return + if getattr(self,'_tangent1',None) is None: + if self.faceType is 'x': t = np.r_[0,0,1.] + elif self.faceType is 'y': t = np.r_[0,0,1.] + elif self.faceType is 'z': t = np.r_[0,1.,0] + self._tangent1 = t + return self._tangent1 + + @property + def normal(self): + if self.faceType is 'x': n = np.r_[1.,0,0] + elif self.faceType is 'y': n = np.r_[0,1.,0] + elif self.faceType is 'z': n = np.r_[0,0,1.] + return n[:self.dim] + + @property + def index(self): + if self.isleaf: return [self.num] + l = [face.index for face in self.children.flatten(order='F')] + # Flatten the list + # e.g. + # [[1,3],[4]] --> [1, 3, 4] + return [item for sublist in l for item in sublist] + + @property + def area(self): + """area of the face""" + return self.sz.prod() + + @property + def length(self): + if self.dim == 3: raise Exception('face.length is not defined for 2D face') + return np.sqrt(((self.node1.x0 - self.node0.x0)**2).sum()) + + def refine(self): + if not self.isleaf: return + self.mesh.isNumbered = False + if self.dim == 2: + self._refine2D() + elif self.dim == 3: + self._refine3D() + + def _refine2D(self): + self.children = np.empty(2,dtype=TreeFace) + # Create refined x0's + x0r_0 = self.x0 + x0r_1 = self.x0+0.5*self.tangent0*self.sz + self.children[0] = TreeFace(self.mesh, x0=x0r_0, faceType=self.faceType, sz=0.5*self.sz, depth=self.depth+1, node0=self.node0) + self.children[1] = TreeFace(self.mesh, x0=x0r_1, faceType=self.faceType, sz=0.5*self.sz, depth=self.depth+1, node0=self.children[0].node1, node1=self.node1) + self.mesh.faces.remove(self) + if self.faceType is 'x': + self.mesh.facesX.remove(self) + elif self.faceType is 'y': + self.mesh.facesY.remove(self) + + def _refine3D(self): + # + # 2_______________3 _______________ + # | e1--> | | | | + # ^ | | ^ | (0,1) | (1,1) | + # | | | | | | | + # | | x | | ---> |-------+-------| + # e2 | | e3 | | | + # | | | (0,0) | (1,0) | + # |_______________| |_______|_______| + # 0 e0--> 1 + + + order = [{'c':[0,0], + 'e0': ('p', 'e0', [0]), 'e1': 'new' , + 'e2': ('p', 'e2', [0]), 'e3': 'new' }, + {'c':[1,0], + 'e0': ('p', 'e0', [1]), 'e1': 'new' , + 'e2': ('c', 'e3', [0,0]), 'e3': ('p', 'e3', [0])}, + {'c':[0,1], + 'e0': ('c', 'e1', [0,0]), 'e1': ('p', 'e1', [0]), + 'e2': ('p', 'e2', [1]), 'e3': 'new' }, + {'c':[1,1], + 'e0': ('c', 'e1', [1,0]), 'e1': ('p', 'e1', [1]), + 'e2': ('c', 'e3', [0,1]), 'e3': ('p', 'e3', [1])}] + + def getEdge(pointer): + if pointer is 'new': return + if pointer[0] == 'p': + return getattr(self, 'edg' + pointer[1]).children[pointer[2][0]] + if pointer[0] == 'c': + f = self.children[pointer[2][0],pointer[2][1]] + return getattr(f, 'edg' + pointer[1]) + + self.children = np.empty((2,2), dtype=TreeFace) + + for edge in [self.edge0, self.edge1, self.edge2, self.edge3]: + edge.refine() + + for O in order: + i, j = O['c'] + x0r = self.x0 + 0.5*i*self.tangent0*self.sz[0] + 0.5*j*self.tangent1*self.sz[1] + e0, e1, e2, e3 = getEdge(O['e0']), getEdge(O['e1']), getEdge(O['e2']), getEdge(O['e3']) + self.children[i,j] = TreeFace(self.mesh, x0=x0r, faceType=self.faceType, depth=self.depth+1, sz=0.5*self.sz, edge0=e0, edge1=e1, edge2=e2, edge3=e3) + + self.mesh.faces.remove(self) + if self.faceType is 'x': + self.mesh.facesX.remove(self) + elif self.faceType is 'y': + self.mesh.facesY.remove(self) + elif self.faceType is 'z': + self.mesh.facesZ.remove(self) + + def plotGrid(self, ax, text=True): + if not self.isleaf: return + if self.dim == 2: + line = np.c_[self.node0.x0, self.node1.x0].T + ax.plot(line[:,0], line[:,1],'r-') + if text: ax.text(self.center[0], self.center[1],self.num) + elif self.dim == 3: + if text: ax.text(self.center[0], self.center[1], self.center[2], self.num) + + @property + def center(self): + if self.dim == 2: + return self.x0 + 0.5*self.tangent0*self.sz[0] + elif self.dim == 3: + return self.x0 + 0.5*self.tangent0*self.sz[0] + 0.5*self.tangent1*self.sz[1] + + +class TreeCell(object): + + __slots__ = ['mesh', 'children', 'depth', 'num', 'sz', + 'node0', 'node1', 'node2', 'node3', + 'node4', 'node5', 'node6', 'node7', + 'fXm', 'fXp', 'fYm', 'fYp', 'fZm', 'fZp', + 'eX0','eX1','eX2','eX3', + 'eY0','eY1','eY2','eY3', + 'eZ0','eZ1','eZ2','eZ3'] + + def __init__(self, mesh, x0=[0,0], depth=0, sz=[1,1], + fXm=None, fXp=None, + fYm=None, fYp=None, + fZm=None, fZp=None): + + self.mesh = mesh + self.depth = depth + + self.sz = np.array(sz, dtype=float) + if self.dim == 2: + # + # 2___________3 + # | fYp | + # | | + # fXm| x |fXp y + # | | ^ + # |___________| |___> x + # 0 fYm 1 + # + n0 = getattr(fXm, 'node0', None) or getattr(fYm, 'node0', None) + n1 = getattr(fXp, 'node0', None) or getattr(fYm, 'node1', None) + n2 = getattr(fXm, 'node1', None) or getattr(fYp, 'node0', None) + n3 = getattr(fXp, 'node1', None) or getattr(fYp, 'node1', None) + + self.fXm = fXm if isinstance(fXm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] ], faceType='x', sz=np.r_[sz[1]], depth=depth, node0=n0, node1=n2) + n0, n2 = self.fXm.node0, self.fXm.node1 + + self.fXp = fXp if isinstance(fXp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0]+sz[0], x0[1] ], faceType='x', sz=np.r_[sz[1]], depth=depth, node0=n1, node1=n3) + n1, n3 = self.fXp.node0, self.fXp.node1 + + self.fYm = fYm if isinstance(fYm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] ], faceType='y', sz=np.r_[sz[0]], depth=depth, node0=n0, node1=n1) + n0, n1 = self.fYm.node0, self.fYm.node1 + + self.fYp = fYp if isinstance(fYp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1]+sz[1]], faceType='y', sz=np.r_[sz[0]], depth=depth, node0=n2, node1=n3) + n2, n3 = self.fYp.node0, self.fYp.node1 + + self.node0, self.node1, self.node2, self.node3 = n0, n1, n2, n3 + + elif self.dim == 3: + # fZp + # | + # 6 ------eX3------ 7 + # /| | / | + # /eZ2 . / eZ3 + # eY2 | fYp eY3 | + # / | / fXp| + # 4 ------eX2----- 5 | + # |fXm 2 -----eX1--|---- 3 z + # eZ0 / | eY1 ^ y + # | eY0 . fYm eZ1 / | / + # | / | | / | / + # 0 ------eX0------1 o----> x + # | + # fZm + # + # + # fX fY fZ + # 2___________3 2___________3 2___________3 + # | e1 | | e1 | | e1 | + # | | | | | | + # e2 | x | e3 z e2 | x | e3 z e2 | x | e3 y + # | | ^ | | ^ | | ^ + # |___________| |___> y |___________| |___> x |___________| |___> x + # 0 e0 1 0 e0 1 0 e0 1 + # + # Mapping Nodes: numOnFace > numOnCell + # + # fXm 0>0, 1>2, 2>4, 3>6 fYm 0>0, 1>1, 2>4, 3>5 fZm 0>0, 1>1, 2>2, 3>3 + # fXp 0>1, 1>3, 2>5, 3>7 fYp 0>2, 1>3, 2>6, 3>7 fZp 0>4, 1>5, 2>6, 3>7 + + def getEdge(face, key): + if face is None: return + return getattr(face, key) + + E = {} + eX0 = getEdge(fYm, 'edge0') or getEdge(fZm, 'edge0') + eX1 = getEdge(fYp, 'edge0') or getEdge(fZm, 'edge1') + eX2 = getEdge(fYm, 'edge1') or getEdge(fZp, 'edge0') + eX3 = getEdge(fYp, 'edge1') or getEdge(fZp, 'edge1') + + eY0 = getEdge(fXm, 'edge0') or getEdge(fZm, 'edge2') + eY1 = getEdge(fXp, 'edge0') or getEdge(fZm, 'edge3') + eY2 = getEdge(fXm, 'edge1') or getEdge(fZp, 'edge2') + eY3 = getEdge(fXp, 'edge1') or getEdge(fZp, 'edge3') + + eZ0 = getEdge(fXm, 'edge2') or getEdge(fYm, 'edge2') + eZ1 = getEdge(fXp, 'edge2') or getEdge(fYm, 'edge3') + eZ2 = getEdge(fXm, 'edge3') or getEdge(fYp, 'edge2') + eZ3 = getEdge(fXp, 'edge3') or getEdge(fYp, 'edge3') + + + self.fXm = fXm if isinstance(fXm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='x', sz=np.r_[sz[1], sz[2]], depth=depth, edge0=eY0, edge1=eY2, edge2=eZ0, edge3=eZ2) + eY0, eY2, eZ0, eZ2 = self.fXm.edge0, self.fXm.edge1, self.fXm.edge2, self.fXm.edge3 + + self.fXp = fXp if isinstance(fXp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0]+sz[0], x0[1] , x0[2] ], faceType='x', sz=np.r_[sz[1], sz[2]], depth=depth, edge0=eY1, edge1=eY3, edge2=eZ1, edge3=eZ3) + eY1, eY3, eZ1, eZ3 = self.fXp.edge0, self.fXp.edge1, self.fXp.edge2, self.fXp.edge3 + + self.fYm = fYm if isinstance(fYm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='y', sz=np.r_[sz[0], sz[2]], depth=depth, edge0=eX0, edge1=eX2, edge2=eZ0, edge3=eZ1) + eX0, eX2, eZ0, eZ1 = self.fYm.edge0, self.fYm.edge1, self.fYm.edge2, self.fYm.edge3 + + self.fYp = fYp if isinstance(fYp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1]+sz[1], x0[2] ], faceType='y', sz=np.r_[sz[0], sz[2]], depth=depth, edge0=eX1, edge1=eX3, edge2=eZ2, edge3=eZ3) + eX1, eX3, eZ2, eZ3 = self.fYp.edge0, self.fYp.edge1, self.fYp.edge2, self.fYp.edge3 + + self.fZm = fZm if isinstance(fZm, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2] ], faceType='z', sz=np.r_[sz[0], sz[1]], depth=depth, edge0=eX0, edge1=eX1, edge2=eY0, edge3=eY1) + eX0, eX1, eY0, eY1 = self.fZm.edge0, self.fZm.edge1, self.fZm.edge2, self.fZm.edge3 + + self.fZp = fZp if isinstance(fZp, TreeFace) else TreeFace(mesh, x0=np.r_[x0[0] , x0[1] , x0[2]+sz[2]], faceType='z', sz=np.r_[sz[0], sz[1]], depth=depth, edge0=eX2, edge1=eX3, edge2=eY2, edge3=eY3) + eX2, eX3, eY2, eY3 = self.fZp.edge0, self.fZp.edge1, self.fZp.edge2, self.fZp.edge3 + + self.eX0, self.eX1, self.eX2, self.eX3, self.eY0, self.eY1, self.eY2, self.eY3, self.eZ0, self.eZ1, self.eZ2, self.eZ3 = eX0, eX1, eX2, eX3, eY0, eY1, eY2, eY3, eZ0, eZ1, eZ2, eZ3 + self.node0, self.node1, self.node2, self.node3, self.node4, self.node5, self.node6, self.node7 = self.fZm.node0, self.fZm.node1, self.fZm.node2, self.fZm.node3, self.fZp.node0, self.fZp.node1, self.fZp.node2, self.fZp.node3 + + mesh.cells.add(self) + + @property + def x0(self): return self.node0.x0 + + @property + def center(self): return self.x0 + 0.5*self.sz + + @property + def dim(self): return self.mesh.dim + + @property + def faceDict(self): + d = {"fXm":self.fXm, "fXp":self.fXp, "fYm":self.fYm, "fYp":self.fYp} + if self.dim == 3: + d["fZm"] = self.fZm + d["fZp"] = self.fZp + return d + + @property + def edgeDict(self): + if self.dim == 2: return None + return {'eX0': self.eX0, 'eX1': self.eX1, 'eX2': self.eX2, 'eX3': self.eX3, 'eY0': self.eY0, 'eY1': self.eY1, 'eY2': self.eY2, 'eY3': self.eY3, 'eZ0': self.eZ0, 'eZ1': self.eZ1, 'eZ2': self.eZ2, 'eZ3': self.eZ3} + + @property + def faceList(self): + l = [self.fXm, self.fXp, self.fYm, self.fYp] + if self.dim == 3: + l += [self.fZm, self.fZp] + return l + + @property + def edgeList(self): + if self.dim == 2: return None + return [self.eX0, self.eX1, self.eX2, self.eX3, self.eY0, self.eY1, self.eY2, self.eY3, self.eZ0, self.eZ1, self.eZ2, self.eZ3] + + @property + def isleaf(self): return getattr(self, 'children', None) is None + + def refine(self, function=None): + if not self.isleaf and function is None: return + + if function is not None: + do = function(self.center) > self.depth + if not do: return + + if self.dim == 2: + self._refine2D() + elif self.dim == 3: + self._refine3D() + + # pass the refine function to the children + if function is not None: + for child in self.children.flatten(): + child.refine(function) + + def _refine2D(self): + + self.mesh.isNumbered = False + + self.children = np.empty((2,2), dtype=TreeCell) + x0, sz = self.x0, self.sz + + for face in self.faceList: + face.refine() + + order = [{'c':[0,0], + 'fXm': ('p', 'fXm', [0]), 'fXp': 'new' , + 'fYm': ('p', 'fYm', [0]), 'fYp': 'new' }, + {'c':[1,0], + 'fXm': ('c', 'fXp', [0,0]), 'fXp': ('p', 'fXp', [0]), + 'fYm': ('p', 'fYm', [1]), 'fYp': 'new' }, + {'c':[0,1], + 'fXm': ('p', 'fXm', [1]), 'fXp': 'new' , + 'fYm': ('c', 'fYp', [0,0]), 'fYp': ('p', 'fYp', [0])}, + {'c':[1,1], + 'fXm': ('c', 'fXp', [0,1]), 'fXp': ('p', 'fXp', [1]), + 'fYm': ('c', 'fYp', [1,0]), 'fYp': ('p', 'fYp', [1])}] + + def getFace(pointer): + if pointer is 'new': return None + if pointer[0] == 'p': + return self.faceDict[pointer[1]].children[pointer[2][0],] + if pointer[0] == 'c': + return self.children[pointer[2][0],pointer[2][1]].faceDict[pointer[1]] + + for O in order: + i, j = O['c'] + x0r = np.r_[x0[0] + 0.5*i*sz[0], x0[1] + 0.5*j*sz[1]] + fXm, fXp, fYm, fYp = getFace(O['fXm']), getFace(O['fXp']), getFace(O['fYm']), getFace(O['fYp']) + self.children[i,j] = TreeCell(self.mesh, x0=x0r, depth=self.depth+1, sz=0.5*sz, fXm=fXm, fXp=fXp, fYm=fYm, fYp=fYp) + + self.mesh.cells.remove(self) + + + def _refine3D(self): + # .----------------.----------------. + # /| /| /| + # / | / | / | + # / | 011 / | 111 / | + # / | / | / | + # .----------------.----+-----------. | + # /| . ---------/|----.----------/|----. + # / | /| / | /| / | /| + # / | / | 001 / | / | 101 / | / | + # / | / | / | / | / | / | + # . -------------- .----------------. |/ | + # | . ---+------|----.----+------|----. | + # | /| .______|___/|____.______|___/|____. + # | / | / 010 | / | / 110| / | / + # | / | / | / | / | / | / + # . ---+---------- . ---+---------- . | / + # | |/ | |/ | |/ z + # | . ----------|----.-----------|----. ^ y + # | / 000 | / 100 | / | / + # | / | / | / | / + # | / | / | / o----> x + # . -------------- . -------------- . + # + # + # Face Refinement: + # + # 2_______________3 _______________ + # | | | | | + # ^ | | | (0,1) | (1,1) | + # | | | | | | + # | | x | ---> |-------+-------| + # t1 | | | | | + # | | | (0,0) | (1,0) | + # |_______________| |_______|_______| + # 0 t0--> 1 + + + order = [{'c':[0,0,0], + 'fXm': ('p', 'fXm', [0,0]), 'fXp': 'new' , + 'fYm': ('p', 'fYm', [0,0]), 'fYp': 'new' , + 'fZm': ('p', 'fZm', [0,0]), 'fZp': 'new' ,}, + {'c':[1,0,0], + 'fXm': ('c', 'fXp', [0,0,0]), 'fXp': ('p', 'fXp', [0,0]), + 'fYm': ('p', 'fYm', [1,0]), 'fYp': 'new' , + 'fZm': ('p', 'fZm', [1,0]), 'fZp': 'new' }, + {'c':[0,1,0], + 'fXm': ('p', 'fXm', [1,0]), 'fXp': 'new' , + 'fYm': ('c', 'fYp', [0,0,0]), 'fYp': ('p', 'fYp', [0,0]), + 'fZm': ('p', 'fZm', [0,1]), 'fZp': 'new' }, + {'c':[1,1,0], + 'fXm': ('c', 'fXp', [0,1,0]), 'fXp': ('p', 'fXp', [1,0]), + 'fYm': ('c', 'fYp', [1,0,0]), 'fYp': ('p', 'fYp', [1,0]), + 'fZm': ('p', 'fZm', [1,1]), 'fZp': 'new' }, + {'c':[0,0,1], + 'fXm': ('p', 'fXm', [0,1]), 'fXp': 'new' , + 'fYm': ('p', 'fYm', [0,1]), 'fYp': 'new' , + 'fZm': ('c', 'fZp', [0,0,0]), 'fZp': ('p', 'fZp', [0,0])}, + {'c':[1,0,1], + 'fXm': ('c', 'fXp', [0,0,1]), 'fXp': ('p', 'fXp', [0,1]), + 'fYm': ('p', 'fYm', [1,1]), 'fYp': 'new' , + 'fZm': ('c', 'fZp', [1,0,0]), 'fZp': ('p', 'fZp', [1,0])}, + {'c':[0,1,1], + 'fXm': ('p', 'fXm', [1,1]), 'fXp': 'new' , + 'fYm': ('c', 'fYp', [0,0,1]), 'fYp': ('p', 'fYp', [0,1]), + 'fZm': ('c', 'fZp', [0,1,0]), 'fZp': ('p', 'fZp', [0,1])}, + {'c':[1,1,1], + 'fXm': ('c', 'fXp', [0,1,1]), 'fXp': ('p', 'fXp', [1,1]), + 'fYm': ('c', 'fYp', [1,0,1]), 'fYp': ('p', 'fYp', [1,1]), + 'fZm': ('c', 'fZp', [1,1,0]), 'fZp': ('p', 'fZp', [1,1])}] + + self.mesh.isNumbered = False + + self.children = np.empty((2,2,2), dtype=TreeCell) + x0, sz = self.x0, self.sz + + for face in self.faceList: + face.refine() + + def getFace(pointer): + if pointer is 'new': return None + if pointer[0] == 'p': + return self.faceDict[pointer[1]].children[pointer[2][0],pointer[2][1]] + if pointer[0] == 'c': + return self.children[pointer[2][0],pointer[2][1],pointer[2][2]].faceDict[pointer[1]] + + for O in order: + i, j, k = O['c'] + x0r = np.r_[x0[0] + 0.5*i*sz[0], x0[1] + 0.5*j*sz[1], x0[2] + 0.5*k*sz[2]] + fXm, fXp, fYm, fYp, fZm, fZp = getFace(O['fXm']), getFace(O['fXp']), getFace(O['fYm']), getFace(O['fYp']), getFace(O['fZm']), getFace(O['fZp']) + self.children[i,j,k] = TreeCell(self.mesh, x0=x0r, depth=self.depth+1, sz=0.5*sz, fXm=fXm, fXp=fXp, fYm=fYm, fYp=fYp, fZm=fZm, fZp=fZp) + + self.mesh.cells.remove(self) + + @property + def faceIndex(self): + F = {} + for face in self.faces: + F[face] = self.faces[face].index + return F + + @property + def vol(self): return self.sz.prod() + + def viz(self, ax, color='none', text=False): + if not self.isleaf: return + x0, sz = self.x0, self.sz + ax.add_patch(plt.Rectangle((x0[0], x0[1]), sz[0], sz[1], facecolor=color, edgecolor='k')) + if text: ax.text(self.center[0],self.center[1],self.num) + + def plotGrid(self, ax, text=False): + if not self.isleaf: return + if self.dim == 2: + ax.plot(self.center[0],self.center[1],'b.') + if text: ax.text(self.center[0],self.center[1],self.num) + elif self.dim == 3: + ax.plot([self.center[0]],[self.center[1]],'b.', zs=[self.center[2]]) + if text: ax.text(self.center[0], self.center[1], self.center[2], self.num) + + +class TreeMesh(InnerProducts, BaseMesh): + """TreeMesh""" + + _meshType = 'TREE' + + def __init__(self, h_in, x0=None): + assert type(h_in) is list, 'h_in must be a list' + h = range(len(h_in)) + 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) + 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. + self.h = h + + if x0 is None: + x0 = np.zeros(self.dim) + else: + assert type(x0) in [list, np.ndarray], 'x0 must be a numpy array or a list' + x0 = np.array(x0, dtype=float) + assert len(x0) == self.dim, 'x0 must have the same dimensions as the mesh' + + # TODO: this has a lot of stuff which doesn't work for this style of mesh... + BaseMesh.__init__(self, np.array([x.size for x in h]), x0) + + # set the sets for holding the cells, nodes, faces, and edges + self.cells = set() + self.nodes = set() + self.faces = set() + self.facesX = set() + self.facesY = set() + if self.dim == 3: + self.facesZ = set() + self.edges = set() + self.edgesX = set() + self.edgesY = set() + self.edgesZ = set() + + self.children = np.empty([hi.size for hi in h],dtype=TreeCell) + + if self.dim == 2: + for i in range(h[0].size): + for j in range(h[1].size): + fXm = None if i is 0 else self.children[i-1][j].fXp + fYm = None if j is 0 else self.children[i][j-1].fYp + x0i = (np.r_[x0[0], h[0][:i]]).sum() + x0j = (np.r_[x0[1], h[1][:j]]).sum() + self.children[i][j] = TreeCell(self, x0=[x0i, x0j], depth=0, sz=[h[0][i], h[1][j]], fXm=fXm, fYm=fYm) + + elif self.dim == 3: + for i in range(h[0].size): + for j in range(h[1].size): + for k in range(h[2].size): + fXm = None if i is 0 else self.children[i-1][j][k].fXp + fYm = None if j is 0 else self.children[i][j-1][k].fYp + fZm = None if k is 0 else self.children[i][j][k-1].fZp + x0i = (np.r_[x0[0], h[0][:i]]).sum() + x0j = (np.r_[x0[1], h[1][:j]]).sum() + x0k = (np.r_[x0[2], h[2][:k]]).sum() + self.children[i][j][k] = TreeCell(self, x0=[x0i, x0j, x0k], depth=0, sz=[h[0][i], h[1][j], h[2][k]], fXm=fXm, fYm=fYm, fZm=fZm) + + isNumbered = Utils.dependentProperty('_isNumbered', False, ['_faceDiv'], 'Setting this to False will delete all operators.') + + @property + def branchdepth(self): + return np.max([node.branchdepth for node in self.children.flatten('F')]) + + def refine(self, function): + for cell in self.children.flatten(): + cell.refine(function) + + def number(self): + if self.isNumbered: return + + self.sortedCells = sorted(self.cells,key=SortByX0()) + for i, sC in enumerate(self.sortedCells): sC.num = i + + self.sortedNodes = sorted(self.nodes,key=SortByX0()) + for i, sN in enumerate(self.sortedNodes): sN.num = i + + self.sortedFaceX = sorted(self.facesX,key=SortByX0()) + for i, sFx in enumerate(self.sortedFaceX): sFx.num = i + + self.sortedFaceY = sorted(self.facesY,key=SortByX0()) + for i, sFy in enumerate(self.sortedFaceY): sFy.num = i + self.nFx + + if self.dim == 3: + self.sortedFaceZ = sorted(self.facesZ,key=SortByX0()) + for i, sFz in enumerate(self.sortedFaceZ): sFz.num = i + self.nFx + self.nFy + + self.sortedEdgeX = sorted(self.edgesX,key=SortByX0()) + for i, sEx in enumerate(self.sortedEdgeX): sEx.num = i + + self.sortedEdgeY = sorted(self.edgesY,key=SortByX0()) + for i, sEy in enumerate(self.sortedEdgeY): sEy.num = i + self.nEx + + self.sortedEdgeZ = sorted(self.edgesZ,key=SortByX0()) + for i, sEz in enumerate(self.sortedEdgeZ): sEz.num = i + self.nEx + self.nEy + + self.isNumbered = True + + @property + def dim(self): return len(self.h) + + @property + def nC(self): return len(self.cells) + + @property + def nN(self): return len(self.nodes) + + @property + def nF(self): return len(self.faces) + + @property + def nFx(self): return len(self.facesX) + + @property + def nFy(self): return len(self.facesY) + + @property + def nFz(self): return None if self.dim < 3 else len(self.facesZ) + + @property + def nE(self): + if self.dim == 2: + return len(self.faces) + elif self.dim == 3: + return len(self.edges) + + @property + def nEx(self): + if self.dim == 2: + return len(self.facesY) + elif self.dim == 3: + return len(self.edgesX) + + @property + def nEy(self): + if self.dim == 2: + return len(self.facesX) + elif self.dim == 3: + return len(self.edgesY) + + @property + def nEz(self): return None if self.dim < 3 else len(self.edgesZ) + + def _grid(self, key): + self.number() + sObjs = {'CC':self.sortedCells, + 'N':self.sortedNodes, + 'Fx': self.sortedFaceX, + 'Fy': self.sortedFaceY, + 'Fz': getattr(self,'sortedFaceZ', None), + 'Ex': getattr(self,'sortedEdgeX', self.sortedFaceY), + 'Ey': getattr(self,'sortedEdgeY', self.sortedFaceX), + 'Ez': getattr(self,'sortedEdgeZ', None)}[key] + G = np.empty((len(sObjs),self.dim)) + for ii, obj in enumerate(sObjs): + G[ii,:] = obj.center + return G + + @property + def gridCC(self): + if getattr(self, '_gridCC', None) is None: + self._gridCC = self._grid('CC') + return self._gridCC + + @property + def gridN(self): + if getattr(self, '_gridN', None) is None: + self._gridN = self._grid('N') + return self._gridN + + @property + def gridFx(self): + if getattr(self, '_gridFx', None) is None: + self._gridFx = self._grid('Fx') + return self._gridFx + + @property + def gridFy(self): + if getattr(self, '_gridFy', None) is None: + self._gridFy = self._grid('Fy') + return self._gridFy + + @property + def gridFz(self): + if self.dim == 2: return None + if getattr(self, '_gridFz', None) is None: + self._gridFz = self._grid('Fz') + return self._gridFz + + @property + def gridEx(self): + if self.dim == 2: return self.gridFy + if getattr(self, '_gridEx', None) is None: + self._gridEx = self._grid('Ex') + return self._gridEx + + @property + def gridEy(self): + if self.dim == 2: return self.gridFx + if getattr(self, '_gridEy', None) is None: + self._gridEy = self._grid('Ey') + return self._gridEy + + @property + def gridEz(self): + if self.dim == 2: return None + if getattr(self, '_gridEz', None) is None: + self._gridEz = self._grid('Ez') + return self._gridEz + + @property + def vol(self): + self.number() + return np.array([cell.vol for cell in self.sortedCells]) + + @property + def area(self): + self.number() + faces = self.sortedFaceX + self.sortedFaceY + if self.dim == 3: + faces += self.sortedFaceZ + return np.array([face.area for face in faces], dtype=float) + + @property + def edge(self): + self.number() + if self.dim == 2: + edges = self.sortedFaceY + self.sortedFaceX + elif self.dim == 3: + edges = self.sortedEdgeX + self.sortedEdgeY + self.sortedEdgeZ + return np.array([e.length for e in edges], dtype=float) + + @property + def faceDiv(self): + if getattr(self, '_faceDiv', None) is None: + self.number() + # TODO: Preallocate! + I, J, V = [], [], [] + for cell in self.sortedCells: + faces = cell.faceDict + for face in faces: + j = faces[face].index + I += [cell.num]*len(j) + J += j + V += [-1 if 'm' in face else 1]*len(j) + VOL = self.vol + D = sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + S = self.area + self._faceDiv = Utils.sdiag(1/VOL)*D*Utils.sdiag(S) + return self._faceDiv + + @property + def edgeCurl(self): + """Construct the 3D curl operator.""" + assert self.dim > 2, "Edge Curl only programed for 3D." + + if getattr(self, '_edgeCurl', None) is None: + self.number() + # TODO: Preallocate! + I, J, V = [], [], [] + for face in self.faces: + for ii, edge in enumerate([face.edge0, face.edge1, face.edge2, face.edge3]): + j = edge.index + I += [face.num]*len(j) + J += j + isNeg = [True, False, True, False] + V += [-1 if isNeg[ii] else 1]*len(j) + C = sp.csr_matrix((V,(I,J)), shape=(self.nF, self.nE)) + S = self.area + L = self.edge + self._edgeCurl = Utils.sdiag(1/S)*C*Utils.sdiag(L) + return self._edgeCurl + + @property + def nodalGrad(self): + if getattr(self, '_nodalGrad', None) is None: + self.number() + # TODO: Preallocate! + I, J, V = [], [], [] + # kinda a hack for the 2D gradient + # because edges are not stored + edges = self.faces if self.dim == 2 else self.edges + for edge in edges: + if self.dim == 3: + I += [edge.num, edge.num] + elif self.dim == 2 and edge.faceType == 'x': + I += [edge.num + self.nFy, edge.num + self.nFy] + elif self.dim == 2 and edge.faceType == 'y': + I += [edge.num - self.nFx, edge.num - self.nFx] + J += [edge.node0.num, edge.node1.num] + V += [-1, 1] + G = sp.csr_matrix((V,(I,J)), shape=(self.nE, self.nN)) + L = self.edge + self._nodalGrad = Utils.sdiag(1/L)*G + return self._nodalGrad + + def _getFaceP(self, face0, face1, face2): + I, J, V = [], [], [] + for cell in self.sortedCells: + face = cell.faceDict[face0] + if face.isleaf: + j = face.index + elif self.dim == 2: + j = face.children[0 if 'm' in face1 else 1].index + elif self.dim == 3: + j = face.children[0 if 'm' in face1 else 1, + 0 if 'm' in face2 else 1].index + lenj = len(j) + I += [cell.num]*lenj + J += j + V += [1./lenj]*lenj + return sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nF)) + + def _getEdgeP(self, edge0, edge1, edge2): + I, J, V = [], [], [] + for cell in self.sortedCells: + if self.dim == 2: + e2f = lambda e: ('f' + {'X':'Y','Y':'X'}[e[1]] + + {'0':'m','1':'p'}[e[2]]) + face = cell.faceDict[e2f(edge0)] + if face.isleaf: + j = face.index + else: + j = face.children[0 if 'm' in e2f(edge1) else 1].index + # Need to flip the numbering for edges + if 'X' in edge0: + j = [jj - self.nFx for jj in j] + elif 'Y' in edge0: + j = [jj + self.nFy for jj in j] + elif self.dim == 3: + edge = cell.edgeDict[edge0] + if edge.isleaf: + j = edge.index + else: + mSide = lambda e: {'0':True,'1':True,'2':False,'3':False}[e[2]] + j = edge.children[0 if mSide(edge1) else 1, + 0 if mSide(edge2) else 1].index + lenj = len(j) + I += [cell.num]*lenj + J += j + V += [1./lenj]*lenj + return sp.csr_matrix((V,(I,J)), shape=(self.nC, self.nE)) + + def _getFacePxx(self, xFace, yFace): + self.number() + xP = self._getFaceP(xFace, yFace, None) + yP = self._getFaceP(yFace, xFace, None) + return sp.vstack((xP, yP)) + + def _getEdgePxx(self, xEdge, yEdge): + self.number() + xP = self._getEdgeP(xEdge, yEdge, None) + yP = self._getEdgeP(yEdge, xEdge, None) + return sp.vstack((xP, yP)) + + def _getFacePxxx(self, xFace, yFace, zFace): + self.number() + xP = self._getFaceP(xFace, yFace, zFace) + yP = self._getFaceP(yFace, xFace, zFace) + zP = self._getFaceP(zFace, xFace, yFace) + return sp.vstack((xP, yP, zP)) + + def _getEdgePxxx(self, xEdge, yEdge, zEdge): + self.number() + xP = self._getEdgeP(xEdge, yEdge, zEdge) + yP = self._getEdgeP(yEdge, xEdge, zEdge) + zP = self._getEdgeP(zEdge, xEdge, yEdge) + return sp.vstack((xP, yP, zP)) + + def plotGrid(self, ax=None, text=True, plotC=True, plotF=True, plotE=False, plotEx=False, plotEy=False, plotEz=False, showIt=False): + axOpts = {'projection':'3d'} if self.dim == 3 else {} + if ax is None: ax = plt.subplot(111, **axOpts) + + if plotC: [c.plotGrid(ax, text=text) for c in self.cells] + if plotF: [f.plotGrid(ax, text=text) for f in self.faces] + if plotE and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edges] + if plotEx and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesX] + if plotEy and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesY] + if plotEz and self.dim==3: [e.plotGrid(ax, text=text) for e in self.edgesZ] + + ax.set_xlim((self.x0[0], self.h[0].sum())) + ax.set_ylim((self.x0[1], self.h[1].sum())) + if self.dim == 3: + ax.set_zlim((self.x0[2], self.h[2].sum())) + if showIt: plt.show() + + def plotImage(self, I, ax=None, showIt=True): + if self.dim == 2: + self._plotImage2D(I, ax=ax, showIt=showIt) + elif self.dim == 3: + raise NotImplementedError('3D visualization is not yet implemented.') + + def _plotImage2D(self, I, ax=None, showIt=True): + if ax is None: ax = plt.subplot(111) + jet = cm = plt.get_cmap('jet') + cNorm = colors.Normalize(vmin=I.min(), vmax=I.max()) + scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet) + ax.set_xlim((self.x0[0], self.h[0].sum())) + ax.set_ylim((self.x0[1], self.h[1].sum())) + for ii, node in enumerate(self.sortedCells): + node.viz(ax=ax, color=scalarMap.to_rgba(I[ii])) + scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots + plt.colorbar(scalarMap) + if showIt: plt.show() + + + +if __name__ == '__main__': + M = TreeMesh([np.ones(x) for x in [4,10]]) + + def function(xc): + r = xc - np.r_[2.,6.] + dist = np.sqrt(r.dot(r)) + if dist < 1.0: + return 3 + if dist < 1.5: + return 2 + else: + return 1 + + M.refine(function) + + DIV = M.faceDiv + Mf = M.getFaceInnerProduct() + # plt.subplot(211) + # plt.spy(DIV) + M.plotGrid(ax=plt.subplot(111),text=True,showIt=True) + + q = np.zeros(M.nC) + q[208] = -1.0 + q[291] = 1.0 + b = Solver(-DIV*Mf*DIV.T).solve(q) + plt.figure() + M.plotImage(b) + # plt.gca().invert_yaxis() + print M.vol + plt.show() diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index 3b8e1eef..3da22a01 100644 --- a/SimPEG/Mesh/__init__.py +++ b/SimPEG/Mesh/__init__.py @@ -1,5 +1,6 @@ from Cyl1DMesh import Cyl1DMesh from TensorMesh import TensorMesh +from TreeMesh import TreeMesh from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh from BaseMesh import BaseMesh from TensorView import TensorView diff --git a/SimPEG/Tests/test_TreeMesh.py b/SimPEG/Tests/test_TreeMesh.py new file mode 100644 index 00000000..f8120859 --- /dev/null +++ b/SimPEG/Tests/test_TreeMesh.py @@ -0,0 +1,503 @@ +from SimPEG.Mesh import TensorMesh +from SimPEG.Mesh.TreeMesh import TreeMesh, TreeFace, TreeCell +import numpy as np +import unittest +import matplotlib.pyplot as plt + +class TestOcTreeObjects(unittest.TestCase): + + def setUp(self): + self.M = TreeMesh([2,1,1]) + self.M.number() + + self.Mr = TreeMesh([2,1,1]) + self.Mr.children[0,0,0].refine() + self.Mr.number() + + def q(s): + if s[0] == 'M': + m = self.M + s = s[1:] + else: + m = self.Mr + c = m.sortedCells[int(s[1])] + if len(s) == 2: return c + if s[2] == 'f' and len(s) == 5: return c.faceDict[s[2:]] + if s[2] == 'f': return getattr(c.faceDict[s[2:5]], 'edg' +s[5:]) + if s[2] == 'e': return getattr(c,s[2:]) + if s[2] == 'n': return getattr(c,'node'+s[3:]) + + self.q = q + + def test_counts(self): + self.assertTrue(self.M.nC == 2) + self.assertTrue(self.M.nFx == 3) + self.assertTrue(self.M.nFy == 4) + self.assertTrue(self.M.nFz == 4) + self.assertTrue(self.M.nF == 11) + self.assertTrue(self.M.nEx == 8) + self.assertTrue(self.M.nEy == 6) + self.assertTrue(self.M.nEz == 6) + self.assertTrue(self.M.nE == 20) + self.assertTrue(self.M.nN == 12) + + self.assertTrue(self.Mr.nC == 9) + self.assertTrue(self.Mr.nFx == 13) + self.assertTrue(self.Mr.nFy == 14) + self.assertTrue(self.Mr.nFz == 14) + self.assertTrue(self.Mr.nF == 41) + + + for cell in self.Mr.sortedCells: + for e in cell.edgeDict: + self.assertTrue(cell.edgeDict[e].edgeType==e[1].lower()) + + self.assertTrue(self.Mr.nN == 31) + self.assertTrue(self.Mr.nEx == 22) + self.assertTrue(self.Mr.nEy == 20) + self.assertTrue(self.Mr.nEz == 20) + + def test_sizes(self): + q = self.q + + for key in ['Mc0','Mc1']: + self.assertTrue(q(key).vol == 0.5) + self.assertTrue(q(key+'fXm').area == 1.) + self.assertTrue(q(key+'fXp').area == 1.) + self.assertTrue(q(key+'fYm').area == 0.5) + self.assertTrue(q(key+'fYp').area == 0.5) + self.assertTrue(q(key+'fZm').area == 0.5) + self.assertTrue(q(key+'fZp').area == 0.5) + + def test_pointersM(self): + q = self.q + + self.assertTrue(q('Mc0fXp') is q('Mc1fXm')) + self.assertTrue(q('Mc0fXpe0') is q('Mc1fXme0')) + self.assertTrue(q('Mc0fXpe1') is q('Mc1fXme1')) + self.assertTrue(q('Mc0fXpe2') is q('Mc1fXme2')) + self.assertTrue(q('Mc0fXpe3') is q('Mc1fXme3')) + self.assertTrue(q('Mc0fYp') is not q('c1fYm')) + self.assertTrue(q('Mc0fXm') is not q('c1fXm')) + + # Test connectivity of shared edges + self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe0')) + self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe1')) + self.assertTrue(q('Mc0fZpe3') is q('Mc1fZpe2')) + self.assertTrue(q('Mc0fZpe3') is not q('c1fZpe3')) + + self.assertTrue(q('Mc0fZme3') is not q('c1fZme0')) + self.assertTrue(q('Mc0fZme3') is not q('c1fZme1')) + self.assertTrue(q('Mc0fZme3') is q('Mc1fZme2')) + self.assertTrue(q('Mc0fZme3') is not q('c1fZme3')) + + self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe0')) + self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe1')) + self.assertTrue(q('Mc0fYpe3') is q('Mc1fYpe2')) + self.assertTrue(q('Mc0fYpe3') is not q('c1fYpe3')) + + self.assertTrue(q('Mc0fYme3') is not q('c1fYme0')) + self.assertTrue(q('Mc0fYme3') is not q('c1fYme1')) + self.assertTrue(q('Mc0fYme3') is q('Mc1fYme2')) + self.assertTrue(q('Mc0fYme3') is not q('c1fYme3')) + + self.assertTrue(q('Mc0fZme3') is q('Mc1fXme0')) + self.assertTrue(q('Mc0fZpe3') is q('Mc1fXme1')) + self.assertTrue(q('Mc0fYme3') is q('Mc1fXme2')) + self.assertTrue(q('Mc0fYpe3') is q('Mc1fXme3')) + + self.assertTrue(q('Mc0fZme3') is q('Mc0fXpe0')) + self.assertTrue(q('Mc0fZpe3') is q('Mc0fXpe1')) + self.assertTrue(q('Mc0fYme3') is q('Mc0fXpe2')) + self.assertTrue(q('Mc0fYpe3') is q('Mc0fXpe3')) + + self.assertTrue(q('Mc1fZme2') is q('Mc1fXme0')) + self.assertTrue(q('Mc1fZpe2') is q('Mc1fXme1')) + self.assertTrue(q('Mc1fYme2') is q('Mc1fXme2')) + self.assertTrue(q('Mc1fYpe2') is q('Mc1fXme3')) + + self.assertTrue(q('Mc1fZme2') is q('Mc0fXpe0')) + self.assertTrue(q('Mc1fZpe2') is q('Mc0fXpe1')) + self.assertTrue(q('Mc1fYme2') is q('Mc0fXpe2')) + self.assertTrue(q('Mc1fYpe2') is q('Mc0fXpe3')) + + + def test_nodePointers(self): + q = self.q + c0 = self.Mr.sortedCells[0] + c0n0 = c0.node0 + self.assertTrue(c0n0 is q('c0n0')) + self.assertTrue(np.all(q('c0n0').center == np.r_[0,0,0.])) + self.assertTrue(q('c0n0').num == 0) + self.assertTrue(q('c0n1').num == 1) + self.assertTrue(q('c0n2').num == 4) + self.assertTrue(q('c0n3').num == 5) + self.assertTrue(q('c0n4').num == 11) + self.assertTrue(q('c0n5').num == 12) + self.assertTrue(q('c0n6').num == 14) + self.assertTrue(q('c0n7').num == 15) + + def test_pointersMr(self): + q = self.q + + c0 = self.Mr.sortedCells[0] + c0fXm = c0.fXm + c0eX0 = c0.eX0 + c0fYme0 = c0.fYm.edge0 + self.assertTrue(c0 is q('c0')) + self.assertTrue(c0fXm is q('c0fXm')) + self.assertTrue(c0eX0 is q('c0eX0')) + self.assertTrue(c0fYme0 is q('c0fYme0')) + + self.assertTrue(q('c0').depth == 1) + self.assertTrue(q('c1').depth == 1) + self.assertTrue(q('c2').depth == 0) + + # Make sure we know where the center of the cells are. + self.assertTrue(np.all(q('c0').center == np.r_[0.125,0.25,0.25])) + self.assertTrue(np.all(q('c1').center == np.r_[0.375,0.25,0.25])) + self.assertTrue(np.all(q('c2').center == np.r_[0.75,0.5,0.5])) + self.assertTrue(np.all(q('c3').center == np.r_[0.125,0.75,0.25])) + self.assertTrue(np.all(q('c4').center == np.r_[0.375,0.75,0.25])) + self.assertTrue(np.all(q('c5').center == np.r_[0.125,0.25,0.75])) + self.assertTrue(np.all(q('c6').center == np.r_[0.375,0.25,0.75])) + self.assertTrue(np.all(q('c7').center == np.r_[0.125,0.75,0.75])) + self.assertTrue(np.all(q('c8').center == np.r_[0.375,0.75,0.75])) + + # Test X face connectivity and locations and stuff... + self.assertTrue(np.all(q('c0fXm').center == np.r_[0,0.25,0.25])) + self.assertTrue(np.all(q('c0fXp').center == np.r_[0.25,0.25,0.25])) + self.assertTrue(q('c0fXp') is q('c1fXm')) + self.assertTrue(np.all(q('c1fXp').center == np.r_[0.5,0.25,0.25])) + self.assertTrue(np.all(q('c2fXm').center == np.r_[0.5,0.5,0.5])) + self.assertTrue(q('c2fXm').branchdepth == 1) + self.assertTrue(q('c2fXm').children[0,0] is q('c1fXp')) + self.assertTrue(np.all(q('c3fXm').center == np.r_[0,0.75,0.25])) + self.assertTrue(np.all(q('c3fXp').center == np.r_[0.25,0.75,0.25])) + self.assertTrue(q('c4fXm') is q('c3fXp')) + self.assertTrue(q('c2fXm').children[1,0] is q('c4fXp')) + + #Test some internal stuff (edges held by cell should be same as inside) + for key in ['Mc0', 'Mc1'] + ['c%d'%i for i in range(9)]: + self.assertTrue(q(key+'eX0') is q(key+'fZme0')) + self.assertTrue(q(key+'eX1') is q(key+'fZme1')) + self.assertTrue(q(key+'eX2') is q(key+'fZpe0')) + self.assertTrue(q(key+'eX3') is q(key+'fZpe1')) + + self.assertTrue(q(key+'eX0') is q(key+'fYme0')) + self.assertTrue(q(key+'eX1') is q(key+'fYpe0')) + self.assertTrue(q(key+'eX2') is q(key+'fYme1')) + self.assertTrue(q(key+'eX3') is q(key+'fYpe1')) + + self.assertTrue(q(key+'eY0') is q(key+'fXme0')) + self.assertTrue(q(key+'eY1') is q(key+'fXpe0')) + self.assertTrue(q(key+'eY2') is q(key+'fXme1')) + self.assertTrue(q(key+'eY3') is q(key+'fXpe1')) + + self.assertTrue(q(key+'eY0') is q(key+'fZme2')) + self.assertTrue(q(key+'eY1') is q(key+'fZme3')) + self.assertTrue(q(key+'eY2') is q(key+'fZpe2')) + self.assertTrue(q(key+'eY3') is q(key+'fZpe3')) + + self.assertTrue(q(key+'eZ0') is q(key+'fXme2')) + self.assertTrue(q(key+'eZ1') is q(key+'fXpe2')) + self.assertTrue(q(key+'eZ2') is q(key+'fXme3')) + self.assertTrue(q(key+'eZ3') is q(key+'fXpe3')) + + self.assertTrue(q(key+'eZ0') is q(key+'fYme2')) + self.assertTrue(q(key+'eZ1') is q(key+'fYme3')) + self.assertTrue(q(key+'eZ2') is q(key+'fYpe2')) + self.assertTrue(q(key+'eZ3') is q(key+'fYpe3')) + + #Test some edge stuff + self.assertTrue(np.all(q('c0eX0').center == np.r_[0.125,0,0])) + self.assertTrue(np.all(q('c0eX1').center == np.r_[0.125,0.5,0])) + self.assertTrue(np.all(q('c0eX2').center == np.r_[0.125,0,0.5])) + self.assertTrue(np.all(q('c0eX3').center == np.r_[0.125,0.5,0.5])) + + self.assertTrue(np.all(q('c5eX0').center == np.r_[0.125,0,0.5])) + self.assertTrue(np.all(q('c5eX1').center == np.r_[0.125,0.5,0.5])) + self.assertTrue(q('c5eX0') is q('c0eX2')) + self.assertTrue(q('c5eX1') is q('c0eX3')) + + self.assertTrue(np.all(q('c0eY0').center == np.r_[0,0.25,0])) + self.assertTrue(np.all(q('c0eY1').center == np.r_[0.25,0.25,0])) + self.assertTrue(np.all(q('c0eY2').center == np.r_[0,0.25,0.5])) + self.assertTrue(np.all(q('c0eY3').center == np.r_[0.25,0.25,0.5])) + + self.assertTrue(np.all(q('c1eY0').center == np.r_[0.25,0.25,0])) + self.assertTrue(np.all(q('c1eY2').center == np.r_[0.25,0.25,0.5])) + self.assertTrue(q('c1eY0') is q('c0eY1')) + self.assertTrue(q('c1eY2') is q('c0eY3')) + + + self.assertTrue(np.all(q('c0eZ0').center == np.r_[0,0,0.25])) + self.assertTrue(np.all(q('c0eZ1').center == np.r_[0.25,0,0.25])) + self.assertTrue(np.all(q('c0eZ2').center == np.r_[0,0.5,0.25])) + self.assertTrue(np.all(q('c0eZ3').center == np.r_[0.25,0.5,0.25])) + + self.assertTrue(np.all(q('c3eZ0').center == np.r_[0,0.5,0.25])) + self.assertTrue(np.all(q('c3eZ1').center == np.r_[0.25,0.5,0.25])) + self.assertTrue(q('c3eZ0') is q('c0eZ2')) + self.assertTrue(q('c3eZ1') is q('c0eZ3')) + + + self.assertTrue(q('c0fXp') is q('c1fXm')) + self.assertTrue(q('c0fYp') is not q('c1fYm')) + self.assertTrue(q('c0fXm') is not q('c1fXm')) + + self.assertTrue(q('c1fXp') is q('c2fXm').children[0,0]) + + self.assertTrue(q('c1fYp') is q('c4fYm')) + self.assertTrue(q('c1fZp') is q('c6fZm')) + + self.assertTrue(q('c6fXp') is q('c2fXm').children[0,1]) + + self.assertTrue(q('c4fXp') is q('c2fXm').children[1,0]) + + + def test_gridCC(self): + x = np.r_[0.25,0.75] + y = np.r_[0.5,0.5] + z = np.r_[0.5,0.5] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridCC).flatten()) == 0) + + x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375] + y = np.r_[0.25,0.25,0.5,0.75,0.75,0.25,0.25,0.75,0.75] + z = np.r_[0.25,0.25,0.5,0.25,0.25,0.75,0.75,0.75,0.75] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridCC).flatten()) == 0) + + def test_gridN(self): + x = np.r_[0,0.5,1,0,0.5,1,0,0.5,1,0,0.5,1] + y = np.r_[0,0,0,1,1,1,0,0,0,1,1,1.] + z = np.r_[0,0,0,0,0,0,1,1,1,1,1,1.] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridN).flatten()) == 0) + + x = np.r_[0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,1] + y = np.r_[0,0,0,0,0.5,0.5,0.5,1,1,1,1,0,0,0,0.5,0.5,0.5,1,1,1,0,0,0,0,0.5,0.5,0.5,1,1,1,1] + z = np.r_[0,0,0,0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1,1,1,1] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridN).flatten()) == 0) + + def test_gridFx(self): + x = np.r_[0.0,0.5,1.0] + y = np.r_[0.5,0.5,0.5] + z = np.r_[0.5,0.5,0.5] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFx).flatten()) == 0) + + x = np.r_[0.0,0.25,0.5,1.0,0.0,0.25,0.5,0.0,0.25,0.5,0.0,0.25,0.5] + y = np.r_[0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.25,0.25,0.25,0.75,0.75,0.75] + z = np.r_[0.25,0.25,0.25,0.5,0.25,0.25,0.25,0.75,0.75,0.75,0.75,0.75,0.75] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFx).flatten()) == 0) + + def test_gridFy(self): + x = np.r_[0.25,0.75,0.25,0.75] + y = np.r_[0,0,1.,1.] + z = np.r_[0.5,0.5,0.5,0.5] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFy).flatten()) == 0) + + x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375] + y = np.r_[0,0,0,0.5,0.5,1,1,1,0,0,0.5,0.5,1,1] + z = np.r_[0.25,0.25,0.5,0.25,0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.75,0.75,0.75] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFy).flatten()) == 0) + + def test_gridFz(self): + x = np.r_[0.25,0.75,0.25,0.75] + y = np.r_[0.5,0.5,0.5,0.5] + z = np.r_[0,0,1.,1.] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridFz).flatten()) == 0) + + x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375,0.125,0.375,0.75,0.125,0.375] + y = np.r_[0.25,0.25,0.5,0.75,0.75,0.25,0.25,0.75,0.75,0.25,0.25,0.5,0.75,0.75] + z = np.r_[0,0,0,0,0,0.5,0.5,0.5,0.5,1,1,1,1,1] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridFz).flatten()) == 0) + + + def test_gridEx(self): + x = np.r_[0.25,0.75,0.25,0.75,0.25,0.75,0.25,0.75] + y = np.r_[0,0,1.,1.,0,0,1.,1.] + z = np.r_[0,0,0,0,1.,1.,1.,1.] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEx).flatten()) == 0) + + x = np.r_[0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.125,0.375,0.125,0.375,0.75,0.125,0.375,0.125,0.375,0.75] + y = np.r_[0,0,0,0.5,0.5,1,1,1,0,0,0.5,0.5,1,1,0,0,0,0.5,0.5,1,1,1] + z = np.r_[0,0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEx).flatten()) == 0) + + def test_gridEy(self): + x = np.r_[0,0.5,1,0,0.5,1] + y = np.r_[0.5,0.5,0.5,0.5,0.5,0.5] + z = np.r_[0,0,0,1.,1.,1.] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEy).flatten()) == 0) + + x = np.r_[0,0.25,0.5,1,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5] + y = np.r_[0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.25,0.25,0.25,0.75,0.75,0.75,0.25,0.25,0.25,0.5,0.75,0.75,0.75] + z = np.r_[0,0,0,0,0,0,0,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEy).flatten()) == 0) + + def test_gridEz(self): + x = np.r_[0,0.5,1,0,0.5,1] + y = np.r_[0,0,0,1.,1.,1.] + z = np.r_[0.5,0.5,0.5,0.5,0.5,0.5] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.M.gridEz).flatten()) == 0) + + x = np.r_[0,0.25,0.5,1,0 ,0.25,0.5,0,0.25,0.5,1,0,0.25,0.5,0 ,0.25,0.5,0 ,0.25,0.5] + y = np.r_[0,0 ,0 ,0,0.5,0.5 ,0.5,1,1 ,1 ,1,0,0 ,0 ,0.5,0.5 ,0.5,1 ,1 ,1 ] + z = np.r_[0.25,0.25,0.25,0.5,0.25,0.25,0.25,0.25,0.25,0.25,0.5,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75,0.75] + self.assertTrue(np.linalg.norm((np.c_[x,y,z]-self.Mr.gridEz).flatten()) == 0) + + +class TestQuadTreeObjects(unittest.TestCase): + + def setUp(self): + self.M = TreeMesh([2,1]) + self.Mr = TreeMesh([2,1]) + self.Mr.children[0,0].refine() + self.Mr.number() + # self.Mr.plotGrid(showIt=True) + + def test_pointersM(self): + c0 = self.M.children[0,0] + c0fXm = c0.fXm + c0fXp = c0.fXp + c0fYm = c0.fYm + c0fYp = c0.fYp + + c1 = self.M.children[1,0] + c1fXm = c1.fXm + c1fXp = c1.fXp + c1fYm = c1.fYm + c1fYp = c1.fYp + + self.assertTrue(c0fXp is c1fXm) + self.assertTrue(c0fYp is not c1fYm) + self.assertTrue(c0fXm is not c1fXm) + + self.assertTrue(c0fXm.area == 1) + self.assertTrue(c0fYm.area == 0.5) + + self.assertTrue(c0.node1 is c1.node0) + self.assertTrue(c0.node3 is c1.node2) + self.assertTrue(self.M.nN == 6) + + + def test_pointersMr(self): + c0 = self.Mr.sortedCells[0] + c0fXm = c0.fXm + c0fXp = c0.fXp + c0fYm = c0.fYm + c0fYp = c0.fYp + + c1 = self.Mr.sortedCells[1] + c1fXm = c1.fXm + c1fXp = c1.fXp + c1fYm = c1.fYm + c1fYp = c1.fYp + + c2 = self.Mr.sortedCells[2] + c2fXm = c2.fXm + c2fXp = c2.fXp + c2fYm = c2.fYm + c2fYp = c2.fYp + + c4 = self.Mr.sortedCells[4] + c4fXm = c4.fXm + c4fXp = c4.fXp + c4fYm = c4.fYm + c4fYp = c4.fYp + + self.assertTrue(c0fXp is c1fXm) + self.assertTrue(c1fXp.node0 is c2fXm.node0) + self.assertTrue(c1fXp.node0 is c2fXm.node0) + self.assertTrue(c4fYm is c1fYp) + self.assertTrue(c4fXp.node1 is c2fXm.node1) + self.assertTrue(c4fXp.node0 is c1fYp.node1) + self.assertTrue(c0fXp.node1 is c4fYm.node0) + + self.assertTrue(self.Mr.nN == 11) + + self.assertTrue(np.all(c1fXp.node0.x0 == np.r_[0.5,0])) + self.assertTrue(np.all(c1fYp.node0.x0 == np.r_[0.25,0.5])) + + +class TestQuadTreeMesh(unittest.TestCase): + + def setUp(self): + M = TreeMesh([np.ones(x) for x in [3,2]]) + for ii in range(1): + M.children[ii,ii].refine() + self.M = M + M.number() + # M.plotGrid(showIt=True) + + def test_MeshSizes(self): + self.assertTrue(self.M.nC==9) + self.assertTrue(self.M.nF==25) + self.assertTrue(self.M.nFx==12) + self.assertTrue(self.M.nFy==13) + self.assertTrue(self.M.nE==25) + self.assertTrue(self.M.nEx==13) + self.assertTrue(self.M.nEy==12) + + def test_gridCC(self): + x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.5,1.5,2.5] + y = np.r_[0.25,0.25,0.5,0.5,0.75,0.75,1.5,1.5,1.5] + self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridCC).flatten()) == 0) + + def test_gridN(self): + x = np.r_[0,0.5,1,2,3,0,0.5,1,0,0.5,1,2,3,0,1,2,3] + y = np.r_[0,0,0,0,0,.5,.5,.5,1,1,1,1,1,2,2,2,2] + self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridN).flatten()) == 0) + + def test_gridFx(self): + x = np.r_[0.0,0.5,1.0,2.0,3.0,0.0,0.5,1.0,0.0,1.0,2.0,3.0] + y = np.r_[0.25,0.25,0.25,0.5,0.5,0.75,0.75,0.75,1.5,1.5,1.5,1.5] + self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridFx).flatten()) == 0) + + def test_gridFy(self): + x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.25,0.75,1.5,2.5,0.5,1.5,2.5] + y = np.r_[0,0,0,0,0.5,0.5,1,1,1,1,2,2,2] + self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridFy).flatten()) == 0) + + def test_gridEx(self): + x = np.r_[0.25,0.75,1.5,2.5,0.25,0.75,0.25,0.75,1.5,2.5,0.5,1.5,2.5] + y = np.r_[0,0,0,0,0.5,0.5,1,1,1,1,2,2,2] + self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridEx).flatten()) == 0) + + def test_gridEy(self): + x = np.r_[0.0,0.5,1.0,2.0,3.0,0.0,0.5,1.0,0.0,1.0,2.0,3.0] + y = np.r_[0.25,0.25,0.25,0.5,0.5,0.75,0.75,0.75,1.5,1.5,1.5,1.5] + self.assertTrue(np.linalg.norm((np.c_[x,y]-self.M.gridEy).flatten()) == 0) + + +class SimpleOctreeOperatorTests(unittest.TestCase): + + def setUp(self): + h1 = np.random.rand(5) + h2 = np.random.rand(7) + h3 = np.random.rand(3) + self.tM = TensorMesh([h1,h2,h3]) + self.oM = TreeMesh([h1,h2,h3]) + self.tM2 = TensorMesh([h1,h2]) + self.oM2 = TreeMesh([h1,h2]) + + def test_faceDiv(self): + self.assertTrue((self.tM.faceDiv - self.oM.faceDiv).toarray().sum() == 0) + self.assertTrue((self.tM2.faceDiv - self.oM2.faceDiv).toarray().sum() == 0) + + def test_nodalGrad(self): + self.assertTrue((self.tM.nodalGrad - self.oM.nodalGrad).toarray().sum() == 0) + self.assertTrue((self.tM2.nodalGrad - self.oM2.nodalGrad).toarray().sum() == 0) + + def test_edgeCurl(self): + self.assertTrue((self.tM.edgeCurl - self.oM.edgeCurl).toarray().sum() == 0) + # self.assertTrue((self.tM2.edgeCurl - self.oM2.edgeCurl).toarray().sum() == 0) + + def test_InnerProducts(self): + self.assertTrue((self.tM.getFaceInnerProduct() - self.oM.getFaceInnerProduct()).toarray().sum() == 0) + self.assertTrue((self.tM2.getFaceInnerProduct() - self.oM2.getFaceInnerProduct()).toarray().sum() == 0) + self.assertTrue((self.tM2.getEdgeInnerProduct() - self.oM2.getEdgeInnerProduct()).toarray().sum() == 0) + self.assertTrue((self.tM.getEdgeInnerProduct() - self.oM.getEdgeInnerProduct()).toarray().sum() == 0) + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 8ed65c84..b2d0d56a 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -136,6 +136,7 @@ def callHooks(match, mainFirst=False): def dependentProperty(name, value, children, doc): def fget(self): return getattr(self,name,value) def fset(self, val): + if getattr(self,name,value) == val: return # it is the same! for child in children: if hasattr(self, child): delattr(self, child)