Skip to content

Commit

Permalink
fix: speeding up interface constraints
Browse files Browse the repository at this point in the history 
few other changes removing nonlinear file, formatting/doc changes
  • Loading branch information
Lachlan Grose committed Nov 1, 2021
1 parent 9136c99 commit 3f4a845
Show file tree
Hide file tree
Showing 4 changed files with 126 additions and 160 deletions.
91 changes: 60 additions & 31 deletions LoopStructural/interpolators/finite_difference_interpolator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,39 +198,68 @@ def add_interface_constraints(self, w=1.0): # for now weight all value points t
# get elements for points
points = self.get_interface_constraints()
if points.shape[0] > 1:
node_idx, inside = self.support.position_to_cell_corners(
points[:, :3])
# print(points[inside,:].shape)

gi = np.zeros(self.support.n_nodes)
gi[:] = -1
gi[self.region] = np.arange(0, self.nx)
idc = np.zeros(node_idx.shape)
idc[:] = -1

idc[inside, :] = gi[node_idx[inside, :]]
inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
a = self.support.position_to_dof_coefs(points[inside, :3]).T
# create oversided array for storing constraints
A = np.zeros((a.shape[0]*a.shape[0],a.shape[1]*2))
interface_idc = np.zeros((a.shape[0]*a.shape[0],a.shape[1]*2),dtype=int)
interface_idc[:] = -1
c_i = 0
vertices, c, tetras, inside = self.support.get_element_for_location(points[:,:3])
# calculate volume of tetras
# vecs = vertices[inside, 1:, :] - vertices[inside, 0, None, :]
# vol = np.abs(np.linalg.det(vecs)) / 6
A = c[inside,:]
# A *= vol[:,None]
idc = tetras[inside,:]

for unique_id in np.unique(points[np.logical_and(~np.isnan(points[:,3]),inside),3]):
mask = points[inside,3] == unique_id
ij = np.array(np.meshgrid(np.arange(0,A[mask,:].shape[0]),np.arange(0,A[mask,:].shape[0]))).T.reshape(-1,2)
interface_A = np.hstack([A[mask,:][ij[:,0],:], -A[mask,:][ij[:,1],:] ])
# interface_A = interface_A.reshape((interface_A.shape[0]*interface_A.shape[0],A.shape[1]))
interface_idc = np.hstack([idc[mask,:][ij[:,0],:], idc[mask,:][ij[:,1],:] ])

# now map the index from global to region create array size of mesh
# initialise as np.nan, then map points inside region to 0->nx
gi = np.zeros(self.support.n_nodes).astype(int)
gi[:] = -1

gi[self.region] = np.arange(0, self.nx)
interface_idc = gi[interface_idc]
# interface_idc = np.tile(interface_idc,(interface_idc.shape[0],1)).reshape(interface_A.shape)#flatten()
outside = ~np.any(interface_idc == -1, axis=1)
self.add_constraints_to_least_squares(interface_A[outside,:] * w,
np.zeros(interface_A[outside,:].shape[0]),
interface_idc[outside, :], name='interface_{}'.format(unique_id))

# if points.shape[0] > 1:
# node_idx, inside = self.support.position_to_cell_corners(
# points[:, :3])
# # print(points[inside,:].shape)

# gi = np.zeros(self.support.n_nodes)
# gi[:] = -1
# gi[self.region] = np.arange(0, self.nx)
# idc = np.zeros(node_idx.shape)
# idc[:] = -1

# idc[inside, :] = gi[node_idx[inside, :]]
# inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
# a = self.support.position_to_dof_coefs(points[inside, :3]).T
# # create oversided array for storing constraints
# A = np.zeros((a.shape[0]*a.shape[0],a.shape[1]*2))
# interface_idc = np.zeros((a.shape[0]*a.shape[0],a.shape[1]*2),dtype=int)
# interface_idc[:] = -1
# c_i = 0

for i in np.unique(points[np.logical_and(~np.isnan(points[:,3]),inside),3]):
mask = points[inside,3] == i
for p1 in range(points[inside][mask].shape[0]):
for p2 in range(p1+1,points[inside][mask].shape[0]):
A[c_i,:8] = a[mask][p1,:]
A[c_i,8:] -= a[mask][p2,:]
interface_idc[c_i,:8] = idc[inside,:][mask,:][p1,:]
interface_idc[c_i,8:] = idc[inside,:][mask,:][p2,:]
c_i+=1
outside = ~np.any(interface_idc == -1, axis=1)
# for i in np.unique(points[np.logical_and(~np.isnan(points[:,3]),inside),3]):
# mask = points[inside,3] == i
# for p1 in range(points[inside][mask].shape[0]):
# for p2 in range(p1+1,points[inside][mask].shape[0]):
# A[c_i,:8] = a[mask][p1,:]
# A[c_i,8:] -= a[mask][p2,:]
# interface_idc[c_i,:8] = idc[inside,:][mask,:][p1,:]
# interface_idc[c_i,8:] = idc[inside,:][mask,:][p2,:]
# c_i+=1
# outside = ~np.any(interface_idc == -1, axis=1)

self.add_constraints_to_least_squares(A[outside,:] * w,
np.zeros(A[outside,:].shape[0]),
interface_idc[outside, :], name='interface')
# self.add_constraints_to_least_squares(A[outside,:] * w,
# np.zeros(A[outside,:].shape[0]),
# interface_idc[outside, :], name='interface')

def add_gradient_constraints(self, w=1.):
"""
Expand Down
58 changes: 0 additions & 58 deletions LoopStructural/interpolators/nonlinear_interpolator.py
View file

This file was deleted.

86 changes: 23 additions & 63 deletions LoopStructural/interpolators/piecewiselinear_interpolator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,6 +254,7 @@ def add_norm_constraints(self, w=1.0):
outside = ~np.any(idc == -1, axis=2)
outside = outside[:, 0]
w = points[:, 6]*w
points[inside,3:6]*=vol[inside]
# w /= 3

self.add_constraints_to_least_squares(d_t[outside, :, :] * w[:,None,None],
Expand Down Expand Up @@ -315,72 +316,31 @@ def add_interface_constraints(self, w=1.0): # for now weight all value points t
points = self.get_interface_constraints()
if points.shape[0] > 1:
vertices, c, tetras, inside = self.support.get_element_for_location(points[:,:3])

gi = np.zeros(self.support.n_nodes)
gi[:] = -1
gi[self.region] = np.arange(0, self.nx)
idc = np.zeros(tetras.shape)
idc[:] = -1

idc[inside, :] = gi[tetras[inside, :]]
inside = np.logical_and(~np.any(idc == -1, axis=1), inside)

# calculate volume of tetras
vecs = vertices[inside, 1:, :] - vertices[inside, 0, None, :]
vol = np.abs(np.linalg.det(vecs)) / 6
a = c[inside]
a *= vol[:,None]
A = c[inside]
A *= vol[:,None]
idc = tetras[inside,:]

# a = self.support.position_to_dof_coefs(points[inside, :3]).T
# create oversided array for storing constraints
A = np.zeros((a.shape[0]*a.shape[0],a.shape[1]*2))
interface_idc = np.zeros((a.shape[0]*a.shape[0],a.shape[1]*2),dtype=int)
interface_idc[:] = -1
c_i = 0

for i in np.unique(points[np.logical_and(~np.isnan(points[:,3]),inside),3]):
mask = points[inside,3] == i
for p1 in range(points[inside][mask].shape[0]):
for p2 in range(p1+1,points[inside][mask].shape[0]):
A[c_i,:4] = a[mask][p1,:]
A[c_i,4:] -= a[mask][p2,:]
interface_idc[c_i,:4] = idc[inside,:][mask,:][p1,:]
interface_idc[c_i,4:] = idc[inside,:][mask,:][p2,:]
c_i+=1
outside = ~np.any(interface_idc == -1, axis=1)

self.add_constraints_to_least_squares(A[outside,:] * w,
np.zeros(A[outside,:].shape[0]),
interface_idc[outside, :], name='interface')
# #get elements for points
# points = self.get_interface_constraints()
# if points.shape[0] > 1:
# vertices, c, tetras, inside = self.support.get_tetra_for_location(points[:,:3])
# # calculate volume of tetras
# vecs = vertices[inside, 1:, :] - vertices[inside, 0, None, :]
# vol = np.abs(np.linalg.det(vecs)) / 6
# A = c[inside]
# A *= vol[:,None]
# idc = tetras[inside,:]
# for unique_id in np.unique(points[np.logical_and(~np.isnan(points[:,3]),inside),3]):
# mask = points[inside,3] == unique_id
# ij = np.array(np.meshgrid(np.arange(0,A[mask,:].shape[0]),np.arange(0,A[mask,:].shape[0]))).T.reshape(-1,2)
# interface_A = np.hstack([A[mask,:][ij[:,0],:], -A[mask,:][ij[:,1],:] ])
# # interface_A = interface_A.reshape((interface_A.shape[0]*interface_A.shape[0],A.shape[1]))
# interface_idc = np.hstack([idc[mask,:][ij[:,0],:], idc[mask,:][ij[:,1],:] ])

# # now map the index from global to region create array size of mesh
# # initialise as np.nan, then map points inside region to 0->nx
# gi = np.zeros(self.support.n_nodes).astype(int)
# gi[:] = -1

# gi[self.region] = np.arange(0, self.nx)
# interface_idc = gi[interface_idc]
# # interface_idc = np.tile(interface_idc,(interface_idc.shape[0],1)).reshape(interface_A.shape)#flatten()
# outside = ~np.any(interface_idc == -1, axis=1)
# self.add_constraints_to_least_squares(interface_A[outside,:] * w,
# np.zeros(interface_A[outside,:].shape[0]),
# interface_idc[outside, :], name='interface_{}'.format(unique_id))
for unique_id in np.unique(points[np.logical_and(~np.isnan(points[:,3]),inside),3]):
mask = points[inside,3] == unique_id
ij = np.array(np.meshgrid(np.arange(0,A[mask,:].shape[0]),np.arange(0,A[mask,:].shape[0]))).T.reshape(-1,2)
interface_A = np.hstack([A[mask,:][ij[:,0],:], -A[mask,:][ij[:,1],:] ])
# interface_A = interface_A.reshape((interface_A.shape[0]*interface_A.shape[0],A.shape[1]))
interface_idc = np.hstack([idc[mask,:][ij[:,0],:], idc[mask,:][ij[:,1],:] ])

# now map the index from global to region create array size of mesh
# initialise as np.nan, then map points inside region to 0->nx
gi = np.zeros(self.support.n_nodes).astype(int)
gi[:] = -1

gi[self.region] = np.arange(0, self.nx)
interface_idc = gi[interface_idc]
# interface_idc = np.tile(interface_idc,(interface_idc.shape[0],1)).reshape(interface_A.shape)#flatten()
outside = ~np.any(interface_idc == -1, axis=1)
self.add_constraints_to_least_squares(interface_A[outside,:] * w,
np.zeros(interface_A[outside,:].shape[0]),
interface_idc[outside, :], name='interface_{}'.format(unique_id))

def add_gradient_orthogonal_constraints(self, points, vector, w=1.0,
B=0):
Expand Down
51 changes: 43 additions & 8 deletions LoopStructural/interpolators/structured_grid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,8 @@ def evaluate_value(self, evaluation_points, property_array):
"""
if property_array.shape[0] != self.n_nodes:
logger.error("Property array does not match grid")
raise BaseException
raise ValueError("cannot assign {} vlaues to array of shape {}".format(
property_array.shape[0], self.n_nodes))
idc, inside = self.position_to_cell_corners(evaluation_points)
v = np.zeros(idc.shape)
v[:, :] = np.nan
Expand All @@ -214,9 +215,35 @@ def evaluate_value(self, evaluation_points, property_array):
return np.sum(v, axis=1)

def evaluate_gradient(self, evaluation_points, property_array):
"""Evaluate the gradient at a location given node values
Parameters
----------
evaluation_points : np.array((N,3))
locations
property_array : np.array((self.nx))
value node, has to be the same length as the number of nodes
Returns
-------
np.array((N,3),dtype=float)
gradient of the implicit function at the locations
Raises
------
ValueError
if the array is not the same shape as the number of nodes
Notes
-----
The implicit function gradient is not normalised, to convert to
a unit vector normalise using vector/=np.linalg.norm(vector,axis=1)[:,None]
"""
if property_array.shape[0] != self.n_nodes:
logger.error("Property array does not match grid")
raise BaseException
raise ValueError("cannot assign {} vlaues to array of shape {}".format(
property_array.shape[0], self.n_nodes))

idc, inside = self.position_to_cell_corners(evaluation_points)
T = np.zeros((idc.shape[0], 3, 8))
T[inside, :, :] = self.get_element_gradient_for_location(evaluation_points[inside, :])[1]
Expand All @@ -232,9 +259,17 @@ def evaluate_gradient(self, evaluation_points, property_array):

def get_element_gradient_for_location(self, pos):
"""
Calculates the gradient matrix at location pos
:param pos: numpy array of location Nx3
:return: Nx3x4 matrix
Get the gradient of the element at the locations.
Parameters
----------
pos : np.array((N,3),dtype=float)
locations
Returns
-------
vertices, gradient, element, inside
[description]
"""
# 6_ _ _ _ 8
# /| /|
Expand All @@ -250,7 +285,7 @@ def get_element_gradient_for_location(self, pos):
T = np.zeros((pos.shape[0], 3, 8))
x, y, z = self.position_to_local_coordinates(pos)
vertices, inside = self.position_to_cell_vertices(pos)
elements = self.global_node_indicies(np.array(self.position_to_cell_index(pos)))
elements,inside = self.position_to_cell_corners(pos)
T[:, 0, 0] = (1 - z) * (y- 1) # v000
T[:, 0, 1] = (1 - y) * (1 - z) # (y[:, 3] - pos[:, 1]) / div
T[:, 0, 2] = -y * (1 - z) # (pos[:, 1] - y[:, 0]) / div
Expand Down Expand Up @@ -296,6 +331,6 @@ def get_element_for_location(self, pos):
[description]
"""
vertices, inside = self.position_to_cell_vertices(pos)
elements = self.global_node_indicies(np.array(self.position_to_cell_index(pos)))
elements, inside = self.position_to_cell_corners(pos)
a = self.position_to_dof_coefs(pos)
return vertices, a, elements, inside
return vertices, a.T, elements, inside

0 comments on commit 3f4a845

Please sign in to comment.