From 079acb1153fb129a37f090caaeb14f63b58d9f80 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 26 Jul 2013 16:38:12 -0700 Subject: [PATCH] volume in 2D by calculating faceInfo (in utils for now.) --- SimPEG/LogicallyOrthogonalMesh.py | 61 +++++++++++++++++++---- SimPEG/utils.py | 83 ++++++++++++++++++++++++++++++- 2 files changed, 131 insertions(+), 13 deletions(-) diff --git a/SimPEG/LogicallyOrthogonalMesh.py b/SimPEG/LogicallyOrthogonalMesh.py index 7369237c..5bb5482b 100644 --- a/SimPEG/LogicallyOrthogonalMesh.py +++ b/SimPEG/LogicallyOrthogonalMesh.py @@ -1,7 +1,7 @@ import numpy as np from BaseMesh import BaseMesh from DiffOperators import DiffOperators -from utils import mkvc, ndgrid +from utils import mkvc, ndgrid, volTetra, indexCube, faceInfo class LogicallyOrthogonalMesh(BaseMesh, DiffOperators): # , LOMGrid @@ -52,21 +52,60 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators): # , LOMGrid gridN = property(**gridN()) # --------------- Geometries --------------------- + # + # + # ------------------- 2D ------------------------- + # + # node(i,j) node(i,j+1) + # A -------------- B + # | | + # | cell(i,j) | + # | I | + # | | + # D -------------- C + # node(i+1,j) node(i+1,j+1) + # + # ------------------- 3D ------------------------- + # + # + # node(i,j,k+1) node(i,j+1,k+1) + # E --------------- F + # /| / | + # / | / | + # / | / | + # node(i,j,k) node(i,j+1,k) + # A -------------- B | + # | H ----------|---- G + # | /cell(i,j) | / + # | / I | / + # | / | / + # D -------------- C + # node(i+1,j,k) node(i+1,j+1,k) def vol(): doc = "Construct cell volumes of the 3D model as 1d array." def fget(self): if(self._vol is None): - vh = self.h - # Compute cell volumes - if(self.dim == 1): - self._vol = mkvc(vh[0]) - elif(self.dim == 2): - # Cell sizes in each direction - self._vol = mkvc(np.outer(vh[0], vh[1])) - elif(self.dim == 3): - # Cell sizes in each direction - self._vol = mkvc(np.outer(mkvc(np.outer(vh[0], vh[1])), vh[2])) + if self.dim == 2: + A, B, C, D = indexCube('ABCD', np.array([self.nNx, self.nNy])) + normal, area, length = faceInfo(np.c_[self.gridN, np.zeros((self.nC, 1))], A, B, C, D) + self._vol = area + elif self.dim == 3: + # Each polyhedron can be decomposed into 5 tetrahedrons + # T1 = [A B D E]; % cutted edge + # T2 = [B E F G]; % cutted edge + # T3 = [B D E G]; % mid + # T4 = [B C D G]; % cutted edge + # T5 = [D E G H]; % cutted edge + A, B, C, D, E, F, G, H = indexCube('ABCDEFGH', np.array([self.nNx, self.nNy, self.nNz])) + + v1 = volTetra(self.gridN, A, B, D, E) # cutted edge + v2 = volTetra(self.gridN, B, E, F, G) # cutted edge + v3 = volTetra(self.gridN, B, D, E, G) # mid + v4 = volTetra(self.gridN, B, C, D, G) # cutted edge + v5 = volTetra(self.gridN, D, E, G, H) # cutted edge + + self._vol = v1 + v2 + v3 + v4 + v5 return self._vol return locals() _vol = None diff --git a/SimPEG/utils.py b/SimPEG/utils.py index 8aa7f27b..37441e8e 100644 --- a/SimPEG/utils.py +++ b/SimPEG/utils.py @@ -126,7 +126,7 @@ def volTetra(xyz, A, B, C, D): def indexCube(nodes, nN): """ - Returns the index of nodes on the mesh. + Returns the index of nodes on the mesh. Input: @@ -167,7 +167,7 @@ def indexCube(nodes, nN): @author Rowan Cockett - Last modified on: 2013/03/07 + Last modified on: 2013/07/26 """ assert type(nodes) == str, "Nodes must be a str variable: e.g. 'ABCD'" @@ -201,6 +201,85 @@ def indexCube(nodes, nN): return out + +def faceInfo(xyz, A, B, C, D, average=True): + """ + function [N] = faceInfo(y,A,B,C,D) + + Returns the averaged normal, area, and edge lengths for a given set of faces. + + If average option is FALSE then N is a cell array {nA,nB,nC,nD} + + + Input: + xyz - X,Y,Z vertex vector + A,B,C,D - vert index of the face (counter clockwize) + + Options: + average - [true]/false, toggles returning all normals or the average + + Output: + N - average face normal or {nA,nB,nC,nD} if average = false + area - average face area + edgeLengths - exact edge Lengths, 4 column vector [AB, BC, CD, DA] + + see also testFaceNormal testFaceArea + + @author Rowan Cockett + + Last modified on: 2013/07/26 + + """ + + # compute normal that is pointing away from you. + # + # A -------A-B------- B + # | | + # | | + # D-A (X) B-C + # | | + # | | + # D -------C-D------- C + + AB = xyz[B, :] - xyz[A, :] + BC = xyz[C, :] - xyz[B, :] + CD = xyz[D, :] - xyz[C, :] + DA = xyz[A, :] - xyz[D, :] + + def cross(X, Y): + return np.c_[X[1, :]*Y[2, :] - X[2, :]*Y[1, :], + X[2, :]*Y[0, :] - X[0, :]*Y[2, :], + X[0, :]*Y[1, :] - X[1, :]*Y[0, :]] + + nA = cross(AB, DA) + nB = cross(BC, AB) + nC = cross(CD, BC) + nD = cross(DA, CD) + + length = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5 + normalize = lambda x: x/np.kron(np.ones((1, x.shape[1]), length(x), 1)) + if average: + # average the normals at each vertex. + N = (nA + nB + nC + nD)/4 # this is intrinsically weighted by area + # normalize + N = normalize(N) + else: + N = [normalize(nA), normalize(nB), normalize(nC), normalize(nD)] + + # Area calculation + # + # Approximate by 4 different triangles, and divide by 2. + # Each triangle is one half of the length of the cross product + # + # So also could be viewed as the average parallelogram. + area = (length(nA)+length(nB)+length(nC)+length(nD))/4 + + # simple edge length calculations + edgeLengths = [length(AB), length(BC), length(CD), length(DA)] + + return N, area, edgeLengths + + def getSubArray(A, ind): """subArray""" return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]]