Skip to content

Commit

Permalink
Merge pull request #660 from Verdant-Evolution/minor-refactoring
Browse files Browse the repository at this point in the history 
Refactored two functions
  • Loading branch information
@kevindlewis23
kevindlewis23 authored Jul 2, 2024
2 parents fd52b54 + b8eee96 commit fb05f8c
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 81 deletions.
65 changes: 23 additions & 42 deletions hexrd/gridutil.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,9 +93,9 @@ def cellIndices(edges, points_1d):
def _fill_connectivity(out, m, n, p):
i_con = 0
for k in range(p):
extra = k*(n+1)*(m+1)
for j in range(m):
for i in range(n):
extra = k*(n+1)*(m+1)
out[i_con, 0] = i + j*(n + 1) + 1 + extra
out[i_con, 1] = i + j*(n + 1) + extra
out[i_con, 2] = i + j + n*(j+1) + 1 + extra
Expand All @@ -122,12 +122,11 @@ def cellConnectivity(m, n, p=1, origin='ul'):

if p > 1:
nele = m*n*(p-1)
tmp_con3 = con.reshape(p, m*n, 4)
tmp_con3 = con.reshape((p, m*n, 4))
hex_con = []
for layer in range(p - 1):
hex_con.append(np.hstack([tmp_con3[layer], tmp_con3[layer + 1]]))
con = np.vstack(hex_con)
pass
if origin.lower().strip() == 'll':
con = con[:, ::-1]
return con
Expand All @@ -152,7 +151,6 @@ def cellCentroids(crd, con):
def compute_areas(xy_eval_vtx, conn):
areas = np.empty(len(conn))
for i in range(len(conn)):
c0, c1, c2, c3 = conn[i]
vtx0x, vtx0y = xy_eval_vtx[conn[i, 0]]
vtx1x, vtx1y = xy_eval_vtx[conn[i, 1]]
v0x, v0y = vtx1x-vtx0x, vtx1y-vtx0y
Expand Down Expand Up @@ -201,24 +199,23 @@ def computeArea(polygon):
triv = np.array([[[0, i - 1], [0, i]] for i in range(2, n_vertices)])

area = 0
for i in range(len(triv)):
for [s1, s2] in triv:
tvp = np.diff(
np.hstack([polygon[triv[i][0], :],
polygon[triv[i][1], :]]), axis=0).flatten()
np.hstack([polygon[s1, :],
polygon[s2, :]]), axis=0).flatten()
area += 0.5 * np.cross(tvp[:2], tvp[2:])
return area


def make_tolerance_grid(bin_width, window_width, num_subdivisions,
adjust_window=False, one_sided=False):
if bin_width > window_width:
bin_width = window_width
bin_width = min(bin_width, window_width)
if adjust_window:
window_width = np.ceil(window_width/bin_width)*bin_width
if one_sided:
ndiv = abs(int(window_width/bin_width))
grid = (np.arange(0, 2*ndiv+1) - ndiv)*bin_width
ndiv = 2*ndiv
ndiv *= 2
else:
ndiv = int(num_subdivisions*np.ceil(window_width/float(bin_width)))
grid = np.arange(0, ndiv+1)*window_width/float(ndiv) - 0.5*window_width
Expand All @@ -228,46 +225,33 @@ def make_tolerance_grid(bin_width, window_width, num_subdivisions,
def computeIntersection(line1, line2):
"""
compute intersection of two-dimensional line intersection
Returns the intersection point as an array of length 2.
If the lines are parallel (or equal) the function returns an empty array.
this is an implementation of two lines:
line1 = [ [x0, y0], [x1, y1] ]
line1 = [ [x3, y3], [x4, y4] ]
line1 = [ [x1, y1], [x2, y2] ]
line2 = [ [x3, y3], [x4, y4] ]
<http://en.wikipedia.org/wiki/Line-line_intersection>
"""
intersection = np.zeros(2)

l1 = np.array(line1)
l2 = np.array(line2)
[x1, y1] = line1[0]
[x2, y2] = line1[1]
[x3, y3] = line2[0]
[x4, y4] = line2[1]

det_l1 = det(l1)
det_l2 = det(l2)
denom = (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
if denom == 0:
return []

det_l1_x = det(np.vstack([l1[:, 0], np.ones(2)]).T)
det_l1_y = det(np.vstack([l1[:, 1], np.ones(2)]).T)
subterm1 = x1*y2 - y1*x2
subterm2 = x3*y4 - y3*x4

det_l2_x = det(np.vstack([l2[:, 0], np.ones(2)]).T)
det_l2_y = det(np.vstack([l2[:, 1], np.ones(2)]).T)

denominator = det(
np.vstack([[det_l1_x, det_l1_y], [det_l2_x, det_l2_y]])
)

if denominator == 0:
intersection = []
else:
intersection[0] = det(
np.vstack(
[[det_l1, det_l1_x], [det_l2, det_l2_x]]
)
) / denominator
intersection[1] = det(
np.vstack(
[[det_l1, det_l1_y], [det_l2, det_l2_y]]
)
) / denominator
intersection[0] = (subterm1*(x3-x4) - subterm2*(x1-x2)) / denom
intersection[1] = (subterm1*(y3-y4) - subterm2*(y1-y2)) / denom
return intersection


Expand Down Expand Up @@ -295,6 +279,7 @@ def isinside(point, boundary, ccw=True):

def sutherlandHodgman(subjectPolygon, clipPolygon):
"""
https://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm
"""
subjectPolygon = np.array(subjectPolygon)
clipPolygon = np.array(clipPolygon)
Expand Down Expand Up @@ -332,7 +317,6 @@ def sutherlandHodgman(subjectPolygon, clipPolygon):
outputList.append(
computeIntersection(subjectLineSegment, clipBoundary)
)
pass
outputList.append(curr_subjectVertex)
elif isinside(prev_subjectVertex, clipBoundary):
subjectLineSegment = np.vstack(
Expand All @@ -341,9 +325,6 @@ def sutherlandHodgman(subjectPolygon, clipPolygon):
outputList.append(
computeIntersection(subjectLineSegment, clipBoundary)
)
pass
prev_subjectVertex = curr_subjectVertex
prev_clipVertex = curr_clipVertex
pass
pass
return outputList
52 changes: 13 additions & 39 deletions hexrd/matrixutil.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,42 +262,19 @@ def vecMVCOBMatrix(R):

T = np.zeros((nrot, 6, 6), dtype='float64')

T[:, 0, 0] = R[:, 0, 0]**2
T[:, 0, 1] = R[:, 0, 1]**2
T[:, 0, 2] = R[:, 0, 2]**2
T[:, 0, 3] = sqr2 * R[:, 0, 1] * R[:, 0, 2]
T[:, 0, 4] = sqr2 * R[:, 0, 0] * R[:, 0, 2]
T[:, 0, 5] = sqr2 * R[:, 0, 0] * R[:, 0, 1]
T[:, 1, 0] = R[:, 1, 0]**2
T[:, 1, 1] = R[:, 1, 1]**2
T[:, 1, 2] = R[:, 1, 2]**2
T[:, 1, 3] = sqr2 * R[:, 1, 1] * R[:, 1, 2]
T[:, 1, 4] = sqr2 * R[:, 1, 0] * R[:, 1, 2]
T[:, 1, 5] = sqr2 * R[:, 1, 0] * R[:, 1, 1]
T[:, 2, 0] = R[:, 2, 0]**2
T[:, 2, 1] = R[:, 2, 1]**2
T[:, 2, 2] = R[:, 2, 2]**2
T[:, 2, 3] = sqr2 * R[:, 2, 1] * R[:, 2, 2]
T[:, 2, 4] = sqr2 * R[:, 2, 0] * R[:, 2, 2]
T[:, 2, 5] = sqr2 * R[:, 2, 0] * R[:, 2, 1]
T[:, 3, 0] = sqr2 * R[:, 1, 0] * R[:, 2, 0]
T[:, 3, 1] = sqr2 * R[:, 1, 1] * R[:, 2, 1]
T[:, 3, 2] = sqr2 * R[:, 1, 2] * R[:, 2, 2]
T[:, 3, 3] = R[:, 1, 2] * R[:, 2, 1] + R[:, 1, 1] * R[:, 2, 2]
T[:, 3, 4] = R[:, 1, 2] * R[:, 2, 0] + R[:, 1, 0] * R[:, 2, 2]
T[:, 3, 5] = R[:, 1, 1] * R[:, 2, 0] + R[:, 1, 0] * R[:, 2, 1]
T[:, 4, 0] = sqr2 * R[:, 0, 0] * R[:, 2, 0]
T[:, 4, 1] = sqr2 * R[:, 0, 1] * R[:, 2, 1]
T[:, 4, 2] = sqr2 * R[:, 0, 2] * R[:, 2, 2]
T[:, 4, 3] = R[:, 0, 2] * R[:, 2, 1] + R[:, 0, 1] * R[:, 2, 2]
T[:, 4, 4] = R[:, 0, 2] * R[:, 2, 0] + R[:, 0, 0] * R[:, 2, 2]
T[:, 4, 5] = R[:, 0, 1] * R[:, 2, 0] + R[:, 0, 0] * R[:, 2, 1]
T[:, 5, 0] = sqr2 * R[:, 0, 0] * R[:, 1, 0]
T[:, 5, 1] = sqr2 * R[:, 0, 1] * R[:, 1, 1]
T[:, 5, 2] = sqr2 * R[:, 0, 2] * R[:, 1, 2]
T[:, 5, 3] = R[:, 0, 2] * R[:, 1, 1] + R[:, 0, 1] * R[:, 1, 2]
T[:, 5, 4] = R[:, 0, 0] * R[:, 1, 2] + R[:, 0, 2] * R[:, 1, 0]
T[:, 5, 5] = R[:, 0, 1] * R[:, 1, 0] + R[:, 0, 0] * R[:, 1, 1]
for i in range(3):
# Other two i values
i1, i2 = [k for k in range(3) if k != i]
for j in range(3):
# Other two j values
j1, j2 = [k for k in range(3) if k != j]

T[:, i, j] = R[:, i, j] ** 2
T[:, i, j + 3] = sqr2 * R[:, i, j1] * R[:, i, j2]
T[:, i + 3, j] = sqr2 * R[:, i1, j] * R[:, i2, j]
T[:, i + 3, j + 3] = (
R[:, i1, j1] * R[:, i2, j2] + R[:, i1, j2] * R[:, i2, j1]
)

if nrot == 1:
T = T.squeeze()
Expand Down Expand Up @@ -562,7 +539,6 @@ def uniqueVectors(v, tol=1.0e-12):
indep = np.hstack([True, tmpcmp > tol]) # independent values
rowint = indep.cumsum()
iv[np.ix_([row], tmpord)] = rowint
pass
#
# Dictionary sort from bottom up
#
Expand All @@ -577,8 +553,6 @@ def uniqueVectors(v, tol=1.0e-12):
if any(ivSrt[:, col] != ivSrt[:, col - 1]):
ivInd[nUniq] = col
nUniq += 1
pass
pass

return vSrt[:, ivInd[0:nUniq]]

Expand Down
53 changes: 53 additions & 0 deletions tests/test_matrix_utils.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
import numpy as np

from hexrd import matrixutil as mutil


def test_vec_mv_cob_matrix():
np.random.seed(0)
# Generate some random matrices
R = np.random.rand(20, 3, 3) * 2 - 1

T = np.zeros((len(R), 6, 6), dtype='float64')
sqr2 = np.sqrt(2)
# Hardcoded implementation
T[:, 0, 0] = R[:, 0, 0]**2
T[:, 0, 1] = R[:, 0, 1]**2
T[:, 0, 2] = R[:, 0, 2]**2
T[:, 0, 3] = sqr2 * R[:, 0, 1] * R[:, 0, 2]
T[:, 0, 4] = sqr2 * R[:, 0, 0] * R[:, 0, 2]
T[:, 0, 5] = sqr2 * R[:, 0, 0] * R[:, 0, 1]
T[:, 1, 0] = R[:, 1, 0]**2
T[:, 1, 1] = R[:, 1, 1]**2
T[:, 1, 2] = R[:, 1, 2]**2
T[:, 1, 3] = sqr2 * R[:, 1, 1] * R[:, 1, 2]
T[:, 1, 4] = sqr2 * R[:, 1, 0] * R[:, 1, 2]
T[:, 1, 5] = sqr2 * R[:, 1, 0] * R[:, 1, 1]
T[:, 2, 0] = R[:, 2, 0]**2
T[:, 2, 1] = R[:, 2, 1]**2
T[:, 2, 2] = R[:, 2, 2]**2
T[:, 2, 3] = sqr2 * R[:, 2, 1] * R[:, 2, 2]
T[:, 2, 4] = sqr2 * R[:, 2, 0] * R[:, 2, 2]
T[:, 2, 5] = sqr2 * R[:, 2, 0] * R[:, 2, 1]
T[:, 3, 0] = sqr2 * R[:, 1, 0] * R[:, 2, 0]
T[:, 3, 1] = sqr2 * R[:, 1, 1] * R[:, 2, 1]
T[:, 3, 2] = sqr2 * R[:, 1, 2] * R[:, 2, 2]
T[:, 3, 3] = R[:, 1, 2] * R[:, 2, 1] + R[:, 1, 1] * R[:, 2, 2]
T[:, 3, 4] = R[:, 1, 2] * R[:, 2, 0] + R[:, 1, 0] * R[:, 2, 2]
T[:, 3, 5] = R[:, 1, 1] * R[:, 2, 0] + R[:, 1, 0] * R[:, 2, 1]
T[:, 4, 0] = sqr2 * R[:, 0, 0] * R[:, 2, 0]
T[:, 4, 1] = sqr2 * R[:, 0, 1] * R[:, 2, 1]
T[:, 4, 2] = sqr2 * R[:, 0, 2] * R[:, 2, 2]
T[:, 4, 3] = R[:, 0, 2] * R[:, 2, 1] + R[:, 0, 1] * R[:, 2, 2]
T[:, 4, 4] = R[:, 0, 2] * R[:, 2, 0] + R[:, 0, 0] * R[:, 2, 2]
T[:, 4, 5] = R[:, 0, 1] * R[:, 2, 0] + R[:, 0, 0] * R[:, 2, 1]
T[:, 5, 0] = sqr2 * R[:, 0, 0] * R[:, 1, 0]
T[:, 5, 1] = sqr2 * R[:, 0, 1] * R[:, 1, 1]
T[:, 5, 2] = sqr2 * R[:, 0, 2] * R[:, 1, 2]
T[:, 5, 3] = R[:, 0, 2] * R[:, 1, 1] + R[:, 0, 1] * R[:, 1, 2]
T[:, 5, 4] = R[:, 0, 0] * R[:, 1, 2] + R[:, 0, 2] * R[:, 1, 0]
T[:, 5, 5] = R[:, 0, 1] * R[:, 1, 0] + R[:, 0, 0] * R[:, 1, 1]

T2 = mutil.vecMVCOBMatrix(R)

assert np.allclose(T, T2)

0 comments on commit fb05f8c

Please sign in to comment.