# 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]

@author Rowan Cockett

"""
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")