bluepulsar
427 lines
14 KiB
Python
427 lines
14 KiB
Python
'''
|
|
Copyright (c) 2015, Harsh Bhatia (bhatia4@llnl.gov)
|
|
All rights reserved.
|
|
|
|
Redistribution and use in source and binary forms, with or without
|
|
modification, are permitted provided that the following conditions are met:
|
|
1. Redistributions of source code must retain the above copyright notice, this
|
|
list of conditions and the following disclaimer.
|
|
2. Redistributions in binary form must reproduce the above copyright notice,
|
|
this list of conditions and the following disclaimer in the documentation
|
|
and/or other materials provided with the distribution.
|
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
|
|
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
|
|
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
|
|
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
|
|
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
|
|
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
|
|
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
|
|
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
|
|
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
|
|
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
'''
|
|
|
|
import sys
|
|
import numpy
|
|
from scipy import spatial
|
|
from .timer import Timer
|
|
|
|
class UnstructuredGrid(object):
|
|
|
|
'''Class to support nHHD on unstructured grids (triangular and tetrahedral)'''
|
|
|
|
# --------------------------------------------------------------------------
|
|
# --------------------------------------------------------------------------
|
|
def need_volumes(self, verbose=False):
|
|
|
|
''' Compute volumes/areas for vertices, simplices, and corners'''
|
|
|
|
if self.pvolumes.shape != (0,):
|
|
return self.pvolumes
|
|
|
|
def tri_area(o,a,b):
|
|
return numpy.linalg.norm(numpy.cross(b-o,a-o)) / 2.0
|
|
|
|
def tet_volume(o,a,b,c):
|
|
return numpy.abs(numpy.dot(a-o,numpy.cross(b-o,c-o))) / 6.0
|
|
|
|
if verbose:
|
|
print(' Computing point areas/volumes...', end ='')
|
|
sys.stdout.flush()
|
|
mtimer = Timer()
|
|
|
|
self.svolumes = numpy.zeros(self.nsimplices,)
|
|
self.pvolumes = numpy.zeros(self.nvertices,)
|
|
self.cvolumes = numpy.zeros((self.nsimplices,self.dim+1))
|
|
|
|
# triangulation
|
|
if self.dim == 2:
|
|
|
|
for sidx in xrange(self.nsimplices):
|
|
|
|
simp = self.simplices[sidx]
|
|
verts = self.vertices[simp]
|
|
|
|
self.svolumes[sidx] = tri_area(verts[0], verts[1], verts[2])
|
|
|
|
# edge weights
|
|
e = [ verts[2]-verts[1], verts[0]-verts[2], verts[1]-verts[0] ]
|
|
l2 = [ numpy.dot(e[0],e[0]), numpy.dot(e[1],e[1]), numpy.dot(e[2],e[2]) ]
|
|
ew = [ l2[0]*(l2[1]+l2[2]-l2[0]), l2[1]*(l2[2]+l2[0]-l2[1]),l2[2]*(l2[0]+l2[1]-l2[2]) ]
|
|
|
|
# corner areas
|
|
if (ew[0] <= 0):
|
|
self.cvolumes[sidx,1] = -0.25 * l2[2] * self.svolumes[sidx] / numpy.dot(e[0], e[2])
|
|
self.cvolumes[sidx,2] = -0.25 * l2[1] * self.svolumes[sidx] / numpy.dot(e[0], e[1])
|
|
self.cvolumes[sidx,0] = self.svolumes[sidx] - self.cvolumes[sidx,1] - self.cvolumes[sidx,2]
|
|
elif (ew[1] <= 0):
|
|
self.cvolumes[sidx,2] = -0.25 * l2[0] * self.svolumes[sidx] / numpy.dot(e[1], e[0]);
|
|
self.cvolumes[sidx,0] = -0.25 * l2[2] * self.svolumes[sidx] / numpy.dot(e[1], e[2]);
|
|
self.cvolumes[sidx,1] = self.svolumes[sidx] - self.cvolumes[sidx,2] - self.cvolumes[sidx,0];
|
|
elif (ew[2] <= 0):
|
|
self.cvolumes[sidx,0] = -0.25 * l2[1] * self.svolumes[sidx] / numpy.dot(e[2], e[1]);
|
|
self.cvolumes[sidx,1] = -0.25 * l2[0] * self.svolumes[sidx] / numpy.dot(e[2], e[0]);
|
|
self.cvolumes[sidx,2] = self.svolumes[sidx] - self.cvolumes[sidx,0] - self.cvolumes[sidx,1];
|
|
else:
|
|
ewscale = 0.5 * self.svolumes[sidx] / (ew[0] + ew[1] + ew[2])
|
|
for d in xrange(3):
|
|
self.cvolumes[sidx,d] = ewscale * (ew[(d+1)%3] + ew[(d+2)%3])
|
|
|
|
self.pvolumes[simp[0]] += self.cvolumes[sidx,0]
|
|
self.pvolumes[simp[1]] += self.cvolumes[sidx,1]
|
|
self.pvolumes[simp[2]] += self.cvolumes[sidx,2]
|
|
|
|
# tetrahedralization
|
|
elif self.sdim == 3:
|
|
|
|
raise ValueError('TODO: pvolumes for 3D')
|
|
for sidx in xrange(self.nsimplices):
|
|
|
|
simp = self.simplices[sidx]
|
|
verts = self.vertices[simp]
|
|
|
|
self.svolumes[sidx] = tet_volume(verts[0], verts[1], verts[2], verts[3])
|
|
|
|
for v in simp:
|
|
self.pvolumes[v] += self.svolumes[sidx] / 4.0
|
|
|
|
|
|
if verbose:
|
|
print(' Done!', end ='')
|
|
mtimer.end()
|
|
|
|
return self.pvolumes
|
|
|
|
def need_adjacentfaces(self, verbose=False):
|
|
|
|
'''
|
|
Find adjacent faces for each vertex
|
|
as list of lists
|
|
'''
|
|
|
|
if len(self.adjacent_faces) != 0:
|
|
return self.adjacent_faces
|
|
|
|
if verbose:
|
|
print(' Computing adjacent_faces...', end ='')
|
|
sys.stdout.flush()
|
|
mtimer = Timer()
|
|
|
|
numadjacentfaces = numpy.zeros(self.nvertices, dtype=int)
|
|
for f in self.simplices:
|
|
for i in xrange(3):
|
|
numadjacentfaces[f[i]] += 1
|
|
|
|
# can be optimized further by avoiding "append"?
|
|
self.adjacent_faces = [[] for _ in xrange(self.nvertices)]
|
|
for fidx in xrange(self.nsimplices):
|
|
for i in xrange(3):
|
|
self.adjacent_faces[self.simplices[fidx, i]].append(fidx)
|
|
|
|
if verbose:
|
|
print('Done!', end ='')
|
|
mtimer.end()
|
|
|
|
return self.adjacent_faces
|
|
|
|
def need_acrossedge(self, verbose=False):
|
|
|
|
'''
|
|
Find adjacent faces for each face (across each edge)
|
|
as ndarray of ints: shape (nsimplex, 3)
|
|
-1 denotes a face on the boundary (no face across edge)
|
|
'''
|
|
if self.across_edge.shape != (0,0):
|
|
return self.across_edge
|
|
|
|
self.need_adjacentfaces(verbose)
|
|
|
|
if verbose:
|
|
print(' Computing across_edge...', end ='')
|
|
sys.stdout.flush()
|
|
mtimer = Timer()
|
|
|
|
self.across_edge = -1 * numpy.ones((self.nsimplices, 3), dtype=numpy.int)
|
|
|
|
for fidx in xrange(self.nsimplices):
|
|
for i in xrange(3):
|
|
|
|
if self.across_edge[fidx, i] != -1:
|
|
continue
|
|
|
|
v1 = self.simplices[fidx, (i+1)%3]
|
|
v2 = self.simplices[fidx, (i+2)%3]
|
|
|
|
for other in self.adjacent_faces[v1]:
|
|
|
|
if other == fidx:
|
|
continue
|
|
|
|
if other not in self.adjacent_faces[v2]:
|
|
continue
|
|
|
|
oface = self.simplices[other]
|
|
j = numpy.where(oface == v1)[0]
|
|
j = (j+1)%3
|
|
|
|
if oface[(j+1)%3] != v2:
|
|
continue
|
|
|
|
self.across_edge[fidx, i] = other
|
|
self.across_edge[other, j] = fidx
|
|
|
|
if verbose:
|
|
print(' Done!', end ='')
|
|
mtimer.end()
|
|
return self.across_edge
|
|
|
|
def need_boundary(self, verbose=False):
|
|
|
|
'''
|
|
Find boundary of the triangulation
|
|
as collection of boundary edeges as an array of [face, k]
|
|
'''
|
|
|
|
if self.bedges.shape != (0,0):
|
|
return self.bedges
|
|
|
|
self.need_acrossedge(verbose)
|
|
|
|
if verbose:
|
|
print(' Computing the boundary...', end ='')
|
|
sys.stdout.flush()
|
|
mtimer = Timer()
|
|
|
|
# find all boundary faces and edges
|
|
bfaces = [fidx for fidx in xrange(self.nsimplices) if -1 in self.across_edge[fidx]]
|
|
bedges = []
|
|
|
|
for fidx in bfaces:
|
|
face = self.simplices[fidx]
|
|
nbrs = self.across_edge[fidx]
|
|
bedges.extend( [[fidx, k] for k in xrange(3) if nbrs[k] == -1] )
|
|
|
|
self.bedges = numpy.array(bedges)
|
|
if verbose:
|
|
print(' Done! found', self.bedges.shape[0], 'boundary edges', end ='')
|
|
mtimer.end()
|
|
|
|
return self.bedges
|
|
|
|
|
|
# --------------------------------------------------------------------------
|
|
# --------------------------------------------------------------------------
|
|
def __init__(self, **kwargs):
|
|
|
|
'''
|
|
kwargs:
|
|
vertices: ndarray of shape (nverts, dim) # dim = 2,3
|
|
simplices: ndarray of shape (nfaces, dim+1)
|
|
verbose: verbosity level
|
|
'''
|
|
|
|
args = kwargs.keys()
|
|
|
|
if 'vertices' not in args:
|
|
raise SyntaxError("Mesh object needs vertex data")
|
|
|
|
verbose = 2
|
|
if 'verbose' in args:
|
|
verbose = kwargs['verbose']
|
|
|
|
self.vertices = kwargs['vertices']
|
|
self.dim = self.vertices.shape[1]
|
|
self.nvertices = self.vertices.shape[0]
|
|
|
|
if self.dim != 2 and self.dim != 3:
|
|
raise SyntaxError("Mesh object works for 2D and 3D only")
|
|
|
|
if verbose > 0:
|
|
print(' Initializing', self.dim, 'D mesh with', self.nvertices, 'vertices...', end ='')
|
|
sys.stdout.flush()
|
|
mtimer = Timer()
|
|
|
|
if verbose > 1:
|
|
print('')
|
|
|
|
# create simplices if needed
|
|
if 'simplices' in args:
|
|
self.Delaunay = None
|
|
self.simplices = kwargs['simplices']
|
|
if verbose > 1:
|
|
print(' got', self.simplices.shape[0], 'simplices')
|
|
else:
|
|
if verbose > 1:
|
|
print(' creating Delaunay mesh...')
|
|
sys.stdout.flush()
|
|
|
|
self.Delaunay = spatial.Delaunay(self.vertices)
|
|
self.simplices = self.Delaunay.simplices
|
|
|
|
if verbose > 1:
|
|
print(' Done! created', self.simplices.shape[0], 'simplices')
|
|
sys.stdout.flush()
|
|
|
|
self.nsimplices = self.simplices.shape[0]
|
|
if self.dim != self.simplices.shape[1]-1:
|
|
pdim = self.dim
|
|
sdim = self.simplices.shape[1]-1
|
|
raise SyntaxError("Dimension mismatch! pdim = "+str(pdim)+" and sdim = "+str(sdim)+" do not match!")
|
|
|
|
self.adjacent_faces = []
|
|
self.across_edge = numpy.empty((0,0))
|
|
self.bedges = numpy.empty((0,0))
|
|
self.pvolumes = numpy.empty(0)
|
|
|
|
self.need_volumes(verbose > 1)
|
|
self.need_boundary(verbose > 1)
|
|
#self.need_meshmatrices(verbose > 1)
|
|
|
|
if verbose > 0:
|
|
print(' Done!', end ='')
|
|
mtimer.end()
|
|
|
|
# --------------------------------------------------------------------------
|
|
# --------------------------------------------------------------------------
|
|
def divcurl(self, vfield, verbose=False):
|
|
|
|
if vfield.shape != (self.nsimplices, self.dim):
|
|
print('UnstructuredGrid.divcurl', vfield.shape, self.dim, self.vertices.shape, self.simplices.shape)
|
|
raise ValueError("UnstructuredGrid requires a valid-dimensional vector field")
|
|
|
|
if verbose:
|
|
print(' Computing divcurl...', end ='')
|
|
sys.stdout.flush()
|
|
mtimer = Timer()
|
|
|
|
div = numpy.zeros(self.nvertices)
|
|
curlw = numpy.zeros(self.nvertices)
|
|
|
|
# for each face
|
|
for sidx in xrange(self.nsimplices):
|
|
|
|
simp = self.simplices[sidx]
|
|
|
|
for k in xrange(3):
|
|
|
|
v = simp[k]
|
|
a = simp[ (k+1)%3 ]
|
|
b = simp[ (k+2)%3 ]
|
|
|
|
# normal and tangent vectors
|
|
tvec = self.vertices[b] - self.vertices[a] # counterclockwise
|
|
nvec = numpy.array([-tvec[1], tvec[0]]) # inward
|
|
|
|
dn = numpy.dot(nvec, vfield[sidx])
|
|
tn = numpy.dot(tvec, vfield[sidx])
|
|
|
|
div[v] += dn
|
|
curlw[v] += tn
|
|
|
|
# fix for boundary edges
|
|
for bedge in self.bedges:
|
|
|
|
sidx = bedge[0]
|
|
eidx = bedge[1]
|
|
|
|
a = self.simplices[sidx][(eidx+1)%3]
|
|
b = self.simplices[sidx][(eidx+2)%3]
|
|
|
|
tvec = self.vertices[b] - self.vertices[a] # counterclockwise
|
|
nvec = numpy.array([-tvec[1], tvec[0]]) # inward
|
|
|
|
dn = numpy.dot(nvec, vfield[sidx])
|
|
dt = numpy.dot(tvec, vfield[sidx])
|
|
|
|
div[a] += dn
|
|
div[b] += dn
|
|
curlw[a] += tn
|
|
curlw[b] += tn
|
|
|
|
div = -0.5 * div
|
|
curlw = 0.5 * curlw
|
|
|
|
if verbose:
|
|
print(' Done!', end ='')
|
|
mtimer.end()
|
|
#print('\tdiv =', div.shape, div.min(), div.max())
|
|
#print('\tcurlw =', curlw.shape, curlw.min(), curlw.max())
|
|
|
|
return (div, curlw)
|
|
|
|
def gradient(self, sfield, verbose=False):
|
|
|
|
if sfield.shape[0] != self.nvertices:
|
|
print('UnstructuredGrid.gradient', self.dim, sfield.shape, self.vertices.shape, self.simplices.shape)
|
|
raise ValueError("UnstructuredGrid requires a valid-dimensional vector field")
|
|
|
|
if verbose:
|
|
print(' Computing gradient...', end ='')
|
|
sys.stdout.flush()
|
|
mtimer = Timer()
|
|
|
|
grad = numpy.zeros((self.nsimplices, self.dim))
|
|
|
|
# for 2D
|
|
for sidx in xrange(self.nsimplices):
|
|
|
|
deb = sidx==0
|
|
|
|
simp = self.simplices[sidx]
|
|
f = 0.5 / self.svolumes[sidx]
|
|
|
|
if deb:
|
|
print(sidx, simp, self.svolumes[sidx], f)
|
|
|
|
for k in xrange(3):
|
|
|
|
v = simp[k]
|
|
a = simp[ (k+1)%3 ]
|
|
b = simp[ (k+2)%3 ]
|
|
|
|
# normal and tangent vectors
|
|
tvec = self.vertices[b] - self.vertices[a] # counterclockwise
|
|
nvec = numpy.array([-tvec[1], tvec[0]]) # inward
|
|
|
|
grad[sidx] += f * sfield[v] * nvec
|
|
|
|
if deb:
|
|
print('\t', k, nvec, sfield[v],grad[sidx])
|
|
|
|
if verbose:
|
|
print(' Done!', end ='')
|
|
mtimer.end()
|
|
return grad
|
|
|
|
def rotated_gradient(self, sfield):
|
|
|
|
rgrad = self.gradient(sfield)
|
|
|
|
rgrad[:,[0, 1]] = rgrad[:,[1, 0]]
|
|
rgrad[:,0] *= -1.0
|
|
|
|
return rgrad
|
|
|
|
# --------------------------------------------------------------------------
|
|
# --------------------------------------------------------------------------
|