# -*- coding: utf-8 -*-
"""
General geometry functions associated with unstructured triangular meshes and
functions used within the WaveVal package.
:Dependencies [External]: numpy
:Dependencies [Internal]:
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
# Non-Standard Python Dependencies
import numpy as np
# Local Module Dependencies
# Other Dependencies
# ----------------------------------------------------------------------------
# GLOBAL VARIABLES
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
[docs]class spatialCoverage:
"""
Base class for setting spatial coverage of data in lat/lon ranges.
"""
def __init__(self, lat_min, lat_max, lon_min, lon_max):
self.lat_range = [lat_min, lat_max]
self.lon_range = [lon_min, lon_max]
# ----------------------------------------------------------------------------
# FUNCTION DEFINITIONS
# ----------------------------------------------------------------------------
# ======================== Triangles ==============================
[docs]def getBarycentricWeights(p1, p2, p3, p):
"""
Calculate Barycentric weights for interpolating data from the triangle
node points *p1*, *p2*, *p3* onto the point *p*.
"""
w = (p2[1]-p3[1])*(p1[0]-p3[0]) + (p3[0]-p2[0])*(p1[1]-p3[1])
w1 = ((p2[1]-p3[1])*(p[0]-p3[0]) + (p3[0]-p2[0])*(p[1]-p3[1])) / w
w2 = ((p3[1]-p1[1])*(p[0]-p3[0]) + (p1[0]-p3[0])*(p[1]-p3[1])) / w
w3 = 1.0 - w1 - w2
return [w1, w2, w3]
[docs]def getTriangleArea(p1, p2, p3):
"""
Calculate the area of a triangle defined by the points *p1*, *p2*, *p3*.
"""
area = 0.5*np.abs(p1[0]*(p2[1]-p3[1]) -
p1[1]*(p2[0]-p3[0]) +
(p2[0]*p3[1]-p3[0]*p2[1]))
return area
[docs]def getTriangleCentroid(p1, p2, p3):
"""
Calculate the position of the centroid for a triangle defined by the
points *p1*, *p2*, *p3*.
"""
centroid = [(p1[0]+p2[0]+p3[0])/3., (p1[1]+p2[1]+p3[1])/3.]
return centroid
[docs]def isInsideTriangle(p1, p2, p3, p):
"""
Determine whether the point *p* lies within the triangle defined by
the points *p1*, *p2*, *p3*.
"""
[w1, w2, w3] = getBarycentricWeights(p1, p2, p3, p)
negValues = np.where(np.asarray([w1, w2, w3]) < 0.)[0]
if (len(negValues) == 0):
isInside = True
else:
isInside = False
return isInside
[docs]def interpTriangle(p1, p2, p3, p, v1, v2, v3):
"""
Interpolate the values *v1*, *v2*, *v3* from the triangle nodes
defined by the points *p1*, *p2*, *p3* onto the location *p* using
the Barycentric wighting method.
"""
[w1, w2, w3] = getBarycentricWeights(p1, p2, p3, p)
vp = (w1*v1 + w2*v2 + w3*v3) / (w1 + w2 + w3)
return vp
[docs]def getNearestNode(X, Y, LocX, LocY):
dr = np.sqrt((X-LocX)**2+(Y-LocY)**2)
NNode = np.where(dr == np.min(dr))[0]
return NNode
[docs]def getTriangle(Mesh, LocX, LocY):
X = Mesh.X.copy()
Y = Mesh.Y.copy()
NNode = getNearestNode(X, Y, LocX, LocY)
# elmOff = np.min(Mesh.elems)
# elems = np.where(Mesh.elems[:, :] == NNode[0]+elmOff)[0]
elems = np.where(Mesh.elems[:, :] == NNode[0])[0]
nelems = np.size(elems)
nodes = Mesh.elems[elems, :]
p = [LocX, LocY]
for elem in np.arange(nelems):
p1 = [X[nodes[elem, 0]], Y[nodes[elem, 0]]]
p2 = [X[nodes[elem, 1]], Y[nodes[elem, 1]]]
p3 = [X[nodes[elem, 2]], Y[nodes[elem, 2]]]
if isInsideTriangle(p1, p2, p3, p):
break
else:
continue
if elem == nelems:
NNodes = []
else:
NNodes = nodes[elem, :].copy()
if len(NNode) > 1:
dn = np.diff(NNode)[0]
nl = len(NNode)
NNodes.resize(3*nl)
for lay in np.arange(nl-1):
indx = np.arange(3)+(lay+1)*3
NNodes[indx] = NNodes[0:3]+(lay+1)*dn
return NNodes
[docs]def getWgtsMatrix(NWghts, nstps):
w = np.asarray(NWghts)
w = w[:, np.newaxis]
cols = np.ones([1, nstps])
wgts = w*cols
return wgts
[docs]def getCentroidWeights(p1, p2, p3, Centroid):
CWghts = getBarycentricWeights(p1, p2, p3, Centroid)
return CWghts
[docs]def getNodesWeights(Mesh, LocX, LocY):
NNodes = getTriangle(Mesh, LocX, LocY)
p1 = [Mesh.X[NNodes[0]], Mesh.Y[NNodes[0]]]
p2 = [Mesh.X[NNodes[1]], Mesh.Y[NNodes[1]]]
p3 = [Mesh.X[NNodes[2]], Mesh.Y[NNodes[2]]]
p = [LocX, LocY]
NLocs = [p1, p2, p3]
NWghts = getBarycentricWeights(p1, p2, p3, p)
return NNodes, NWghts, NLocs
# ====================== Scalar Fields ===========================
# ====================== Vector Fields ===========================
[docs]def pol2cart(r, theta):
"""
Coordinate conversion: Polar to Cartesian
:param r: radial distance
:param theta: angle in degrees
:return: x
:return: y
"""
x = r * np.cos(np.deg2rad(theta,dtype='double'),dtype='double')
y = r * np.sin(np.deg2rad(theta,dtype='double'),dtype='double')
return x,y
[docs]def cart2pol(x, y):
"""
Coordinate conversion: Cartesian to Polar
:param x: Cartesian x-coord
:param y: Cartesian y-coord
:return: r
:return: theta (in degrees)
"""
r = np.sqrt(np.power(x,2.0,dtype='double') +
np.power(y,2.0,dtype='double'),dtype='double')
theta = np.rad2deg(np.arctan2(y,x),dtype='double')
return r, theta
[docs]def pol2cmplx(radii, angles, angularUnits=0):
"""
Coordinate conversion: Polar to Complex
:param radii: radial distance
:param angles: angle
:param anglularUnits: 0 = degress [default], 1 = radians
:return: cmplx
"""
if angularUnits == 0:
cmplx = radii * np.exp(1j*angles)
else:
cmplx = radii * np.exp(1j*np.deg2rad(angles))
return cmplx
[docs]def cmplx2pol(x, angularUnits=0):
"""
Coordinate conversion: Complex to Polar
:param x: complex numbers (x + jy)
:param anglularUnits: 0 = degress [default], 1 = radians
:return: radii
:return: angles
"""
if angularUnits == 0:
radii = abs(x)
angles = np.angle(x)
else:
radii = abs(x)
angles = np.rad2deg(np.angle(x))
return radii, angles
[docs]def cart2cmplx(x, y):
"""
Coordinate conversion: Cartesian to Complex
:param x: Cartesian x-coord
:param y: Cartesian y-coord
:return: cmplx
"""
cmplx = x + 1j*y
return cmplx
[docs]def cmplx2cart(c):
"""
Coordinate conversion: Complex to Cartesian
:param c: complex numbers (x + jy)
:return: x
:return: y
"""
x = np.real(c)
y = np.imag(c)
return x, y
[docs]def rotateVectorField(U, V, W, Theta):
thetaRad = np.deg2rad(Theta)
Rz = np.asarray([[np.cos(thetaRad), -np.sin(thetaRad), 0.0],
[np.sin(thetaRad), np.cos(thetaRad), 0.0],
[0.0, 0.0, 1.0]])
orig_shape = U.shape
Uf = np.reshape(U, np.size(U))
Vf = np.reshape(V, np.size(V))
Wf = np.reshape(W, np.size(W))
vel_rot = np.matmul(Rz, [Uf, Vf, Wf])
Ur = np.reshape(vel_rot[0, :], orig_shape)
Vr = np.reshape(vel_rot[1, :], orig_shape)
Wr = np.reshape(vel_rot[2, :], orig_shape)
return Ur, Vr, Wr