Source code for discretize.utils.curvilinear_utils

import numpy as np
from discretize.utils.matrix_utils import mkvc, ndgrid, sub2ind
from discretize.utils.code_utils import deprecate_function
import warnings


[docs]def volume_tetrahedron(xyz, A, B, C, D): """ Returns the volume for tetrahedras volume specified by the indexes A to D. :param numpy.ndarray xyz: X,Y,Z vertex vector :param numpy.ndarray A,B,C,D: vert index of the tetrahedra :rtype: numpy.ndarray :return: V, volume of the tetrahedra Algorithm https://en.wikipedia.org/wiki/Tetrahedron#Volume .. math:: V = {1 \over 3} A h V = {1 \over 6} | ( a - d ) \cdot ( ( b - d ) ( c - d ) ) | """ AD = xyz[A, :] - xyz[D, :] BD = xyz[B, :] - xyz[D, :] CD = xyz[C, :] - xyz[D, :] V = ( (BD[:, 0] * CD[:, 1] - BD[:, 1] * CD[:, 0]) * AD[:, 2] - (BD[:, 0] * CD[:, 2] - BD[:, 2] * CD[:, 0]) * AD[:, 1] + (BD[:, 1] * CD[:, 2] - BD[:, 2] * CD[:, 1]) * AD[:, 0] ) return V / 6
[docs]def index_cube(nodes, grid_size, n=None): """ Returns the index of nodes on the mesh. Input: nodes - string of which nodes to return. e.g. 'ABCD' grid_size - size of the nodal grid n - number of nodes each i,j,k direction: [ni,nj,nk] Output: index - index in the order asked e.g. 'ABCD' --> (A,B,C,D) TWO DIMENSIONS:: node(i,j) node(i,j+1) A -------------- B | | | cell(i,j) | | I | | | D -------------- C node(i+1,j) node(i+1,j+1) THREE DIMENSIONS:: 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) """ if not isinstance(nodes, str): raise TypeError("Nodes must be a str variable: e.g. 'ABCD'") nodes = nodes.upper() try: dim = len(grid_size) if n is None: n = tuple(x - 1 for x in grid_size) except TypeError: return TypeError("grid_size must be iterable") # Make sure that we choose from the possible nodes. possibleNodes = "ABCD" if dim == 2 else "ABCDEFGH" for node in nodes: if node not in possibleNodes: raise ValueError("Nodes must be chosen from: '{0!s}'".format(possibleNodes)) if dim == 2: ij = ndgrid(np.arange(n[0]), np.arange(n[1])) i, j = ij[:, 0], ij[:, 1] elif dim == 3: ijk = ndgrid(np.arange(n[0]), np.arange(n[1]), np.arange(n[2])) i, j, k = ijk[:, 0], ijk[:, 1], ijk[:, 2] else: raise Exception("Only 2 and 3 dimensions supported.") nodeMap = { "A": [0, 0, 0], "B": [0, 1, 0], "C": [1, 1, 0], "D": [1, 0, 0], "E": [0, 0, 1], "F": [0, 1, 1], "G": [1, 1, 1], "H": [1, 0, 1], } out = () for node in nodes: shift = nodeMap[node] if dim == 2: out += (sub2ind(grid_size, np.c_[i + shift[0], j + shift[1]]).flatten(),) elif dim == 3: out += ( sub2ind( grid_size, np.c_[i + shift[0], j + shift[1], k + shift[2]] ).flatten(), ) return out
[docs]def face_info(xyz, A, B, C, D, average=True, normalize_normals=True, **kwargs): """ function [N] = face_info(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 """ if "normalizeNormals" in kwargs: warnings.warn( "The normalizeNormals keyword argument has been deprecated, please use normalize_normals. " "This will be removed in discretize 1.0.0", DeprecationWarning, ) normalize_normals = kwargs["normalizeNormals"] if not isinstance(average, bool): raise TypeError("average must be a boolean") if not isinstance(normalize_normals, bool): raise TypeError("normalize_normals must be a boolean") # 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: np.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2 + x[:, 2] ** 2) normalize = lambda x: x / np.kron(np.ones((1, x.shape[1])), mkvc(length(x), 2)) 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: if normalize_normals: N = [normalize(nA), normalize(nB), normalize(nC), normalize(nD)] else: N = [nA, nB, nC, 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. # # TODO: This does not compute correctly for concave quadrilaterals area = (length(nA) + length(nB) + length(nC) + length(nD)) / 4 return N, area
volTetra = deprecate_function(volume_tetrahedron, "volTetra", removal_version="1.0.0") indexCube = deprecate_function(index_cube, "indexCube", removal_version="1.0.0") faceInfo = deprecate_function(face_info, "faceInfo", removal_version="1.0.0")