# -*- coding: utf-8 -*-
"""
Functions for the manipulation and interpolation of data on unstructured
triangular meshes.
"""
# ----------------------------------------------------------------------------
# IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import numpy as np
# Non-Standard Python Dependencies
# Local Module Dependencies
# Other Dependencies
#--------------------------------------------------------------------------
# GLOBAL CONSTANTS
#--------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# CLASS DEFINITIONS
# ----------------------------------------------------------------------------
#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------
[docs]def getBarycentricWeights(p1,p2,p3,p):
"""
Calculate Barycentric weights for point p within triangle (p1,p2,p3)
Parameters
----------
p1 : list[ float ]
[x, y] coordinates for triangle corner point p1.
p2 : list[ float ]
[x, y] coordinates for triangle corner point p2.
p3 : list[ float ]
[x, y] coordinates for triangle corner point p3.
p : list[ float ]
[x, y] coordinates for point p within triangle (p1, p2, p3).
Returns
-------
[w1, w2 ,w3] : list
barycentric weights for triangle (p1, p2, p3).
"""
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 area of the triangle defined by the points (p1, p2, p3)
Parameters
----------
p1 : list[ float ]
[x, y] coordinates for triangle corner point p1.
p2 : list[ float ]
[x, y] coordinates for triangle corner point p2.
p3 : list[ float ]
[x, y] coordinates for triangle corner point p3.
Returns
-------
area : float
the area of the triangle (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 centroid location for the triangle defined by the points
(p1, p2, p3)
Parameters
----------
p1 : list[ float ]
[x, y] coordinates for triangle corner point p1.
p2 : list[ float ]
[x, y] coordinates for triangle corner point p2.
p3 : list[ float ]
[x, y] coordinates for triangle corner point p3.
Returns
-------
centroid : list[ float ]
centroid location for the triangle (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 if the point p lies with in the triangle defeind by (p1, p2, p3)
Parameters
----------
p1 : list[ float ]
[x, y] coordinates for triangle corner point p1.
p2 : list[ float ]
[x, y] coordinates for triangle corner point p2.
p3 : list[ float ]
[x, y] coordinates for triangle corner point p3.
p : list[ float ]
[x, y] coordinates for point p.
Returns
-------
isInside : bool
True/False flag defining whether p is inside (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):
"""
Use Barycentic weighting to interpolate variable data (v1, v2, v3) on
triangle corner points (p1, p2, p3) onto the location p.
Parameters
----------
p1 : list[ float ]
[x, y] coordinates for triangle corner point p1.
p2 : list[ float ]
[x, y] coordinates for triangle corner point p2.
p3 : list[ float ]
[x, y] coordinates for triangle corner point p3.
p : list[ float ]
[x, y] coordinates for point p.
v1 : float
variable value at point p1.
v2 : floate
variable value at point p2.
v3 : float
variable value at point p3.
Returns
-------
vp : float
interpolated variable value at point p.
"""
[w1,w2,w3] = getBarycentricWeights(p1,p2,p3,p)
vp = (w1*v1 + w2*v2 + w3*v3) / (w1 + w2 + w3)
return vp
[docs]def getElemAroundLoc(meshx,meshy,meshElems,LocX,LocY,nLys=1):
X = np.transpose(meshx.copy())
Y = np.transpose(meshy.copy())
NNode = getNearestNode(X,Y,LocX,LocY,nLys)
nodeElems = np.where(meshElems[:,:] == NNode[0])[0]
nelems = np.size(nodeElems)
nodes = meshElems[nodeElems,:]
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:
nodeElem = None
else:
nodeElem = nodeElems[elem]
return nodeElem
[docs]def getElemNodes(meshx,meshy,elems,LocX,LocY):
elem = getElemAroundLoc(meshx,meshy,elems,LocX,LocY)
if elem is None:
elemNodes = None
else:
elemNodes = elems[elem]
return elem, elemNodes
[docs]def getNearestNode(meshX,meshY,locX,locY,nLys=1):
nPts = len(np.squeeze(meshX))
dr = np.sqrt(np.power((meshX-locX),2)+np.power((meshY-locY),2))
offsets = np.arange(0,nLys)*nPts
nearestNode = np.zeros((nLys),dtype=int)
#nearestNode[:] = np.transpose(np.where(dr == np.min(dr))[0]+offsets)[0]
nn = np.where(dr == np.min(dr))[0]+offsets;
if len(nn) > nLys:
nearestNode[:] = np.transpose(np.where(dr == np.min(dr))[0]+offsets)[0:nLys]
else:
nearestNode[:] = np.transpose(np.where(dr == np.min(dr))[0]+offsets)
return nearestNode
[docs]def getNodesAroundNode(node,edges):
indx, = np.where((edges[:,0] == node) | (edges[:,1] == node))
nedges = indx.size
nodes = np.unique(edges[indx,:].reshape(nedges*2))
return nodes
[docs]def weightsMatrix(nodeWghts,nStps):
w = np.asarray(nodeWghts)
w = w[:,np.newaxis]
cols = np.ones([1,nStps])
wgtsMat = w*cols
return wgtsMat
[docs]def getCentroidWeights(p1,p2,p3):
"""
Calculate Barycentric weights for interpolating data on the centroid of the
triange defined by the points (p1, p2, p3)
Parameters
----------
p1 : list[ float ]
[x, y] coordinates for triangle corner point p1.
p2 : list[ float ]
[x, y] coordinates for triangle corner point p2.
p3 : list[ float ]
[x, y] coordinates for triangle corner point p3.
Returns
-------
centroidWghts : list[ float ]
centroid Barycentric weights for triangle (p1, p2, p3).
"""
centroid = getTriangleCentroid(p1,p2,p3)
centroidWghts = getBarycentricWeights(p1,p2,p3,centroid)
return centroidWghts
[docs]def getNodesWeights(meshx,meshy,elems,LocX,LocY):
"""
Extract the information required to interpolate mesh data onto a specific
location defined by the point (LocX, LocY)
Parameters
----------
meshx : numpy.array, float
mesh nodes x coordinate values.
meshy : numpy.array, float
mesh nodes y coordinate values.
elems : numpy.ndarray, int
mesh elements defninition matrix (nElems, 3).
LocX : float
x coordinate of location of interest.
LocY : flaat
y coordinate of location of interest.
Returns
-------
NNodes : list[ int ]
node indices of triangular element corner points.
NWghts : list[ float ]
Barycentric weighting values for triangular element corner points.
NLocs : list[ list ]
list of the triangular element corner point locations
(p1[x,y], p2[x,y] ,p3[x,y]).
"""
Elem, NNodes = getElemNodes(meshx,meshy,elems,LocX,LocY)
if Elem is None:
NNodes = None
NWghts = None
NLocs = None
else:
p1 = [meshx[NNodes[0]],meshy[NNodes[0]]]
p2 = [meshy[NNodes[1]],meshy[NNodes[1]]]
p3 = [meshy[NNodes[2]],meshy[NNodes[2]]]
p = [LocX,LocY]
NLocs = [p1,p2,p3]
NWghts = getBarycentricWeights(p1,p2,p3,p)
return Elem, NNodes,NWghts,NLocs
[docs]def gradientAtNode(node,nodes,nodesX,nodesY,var):
iref = np.where(nodes == node)
ipts = np.where(nodes != node)
dx = nodesX[ipts]-nodesX[iref]
dy = nodesY[ipts]-nodesY[iref]
A = np.transpose(np.asarray([dx , dy]))
B = var[ipts]-var[iref]
grad, residual, rank, singular = np.linalg.lstsq(A,B,rcond=None)
return grad, residual, rank, singular
[docs]def meshDensity(areas):
mshDen = 1.0/areas
return mshDen
[docs]def gradation2meshDensity(gradRaster):
mshDen = (4.0/np.sqrt(3))/np.square(gradRaster)
return mshDen
[docs]def gradation2numElems(gradRaster,dx,dy):
dA = dx * dy
mshDen = gradation2meshDensity(gradRaster)
numElems = np.int64(np.sum(mshDen*dA))
return numElems