Skip to content

Commit

Permalink
fix: changed to boolean logic for aabb
Browse files Browse the repository at this point in the history
  • Loading branch information
Lachlan Grose committed Oct 11, 2021
1 parent 53492a2 commit f5a5f9b
Show file tree
Hide file tree
Showing 2 changed files with 98 additions and 103 deletions.
70 changes: 8 additions & 62 deletions LoopStructural/interpolators/structured_tetra.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
import logging

import numpy as np
from LoopStructural.interpolators.cython.dsi_helper import cg, constant_norm, fold_cg
from .base_structured_3d_support import BaseStructuredSupport

from LoopStructural.utils import getLogger
Expand Down Expand Up @@ -149,8 +148,7 @@ def get_tetra_for_location(self, pos):
# get cell corners
xi, yi, zi = self.cell_corner_indexes(c_xi, c_yi, c_zi) # global_index_to_node_index(gi)
# convert to node locations
nodes = self.node_indexes_to_position(xi, yi, zi).T

nodes = np.array(self.node_indexes_to_position(xi, yi, zi)).T
vertices[:, :, even_mask, :] = nodes[:, even_mask, :][self.tetra_mask_even, :, :]
vertices[:, :, ~even_mask, :] = nodes[:, ~even_mask, :][self.tetra_mask, :, :]
# changing order to points, tetra, nodes, coord
Expand All @@ -177,7 +175,7 @@ def get_tetra_for_location(self, pos):
c[:, :, 1] = vb / v
c[:, :, 2] = vc / v
c[:, :, 3] = vd / v

print(c)
# if all coords are +ve then point is inside cell
mask = np.all(c > 0, axis=2)

Expand All @@ -198,11 +196,15 @@ def get_tetra_for_location(self, pos):

tetras[even_mask, :, :] = gi[even_mask, :][:, self.tetra_mask_even]
tetras[~even_mask, :, :] = gi[~even_mask, :][:, self.tetra_mask]
inside = np.logical_and(inside,self.inside(pos))
# inside = np.logical_and s(inside,self.inside(pos))
vertices_return = np.zeros((pos.shape[0],4,3))
vertices_return[:] = np.nan
# set all masks not inside to False
mask[~inside,:] = False
print(inside.astype(int)-np.any(mask,axis=1).astype(int))
print(vertices_return[inside,:,:].shape)
print(vertices[mask,:,:].shape)
print(mask)
vertices_return[inside,:,:] = vertices[mask,:,:]#[mask,:,:]#[inside,:,:]
c_return = np.zeros((pos.shape[0],4))
c_return[:] = np.nan
Expand All @@ -212,62 +214,6 @@ def get_tetra_for_location(self, pos):
tetra_return[inside,:] = tetras[mask,:]
return vertices_return, c_return, tetra_return, inside

def get_constant_gradient(self, region, direction=None):
"""
Get the constant gradient for the specified nodes
Parameters
----------
region : np.array(dtype=bool)
mask of nodes to calculate cg for
Returns
-------
"""
"""
Add the constant gradient regularisation to the system
Parameters
----------
w (double) - weighting of the cg parameter
Returns
-------
"""
if direction is not None:
print('using cg direction')
logger.info("Running constant gradient")
elements_gradients = self.get_element_gradients(np.arange(self.ntetra))
if elements_gradients.shape[0] != direction.shape[0]:
logger.error('Cannot add directional CG, vector field is not the correct length')
return
region = region.astype('int64')

neighbours = self.get_neighbours()
elements = self.get_elements()
idc, c, ncons = fold_cg(elements_gradients, direction, neighbours.astype('int64'), elements.astype('int64'), self.nodes)

idc = np.array(idc[:ncons, :])
c = np.array(c[:ncons, :])
B = np.zeros(c.shape[0])
return c, idc, B
if self.cg is None:
logger.info("Running constant gradient")
elements_gradients = self.get_element_gradients(np.arange(self.ntetra))
region = region.astype('int64')

neighbours = self.get_neighbours()
elements = self.get_elements()
idc, c, ncons = cg(elements_gradients, neighbours.astype('int64'), elements.astype('int64'), self.nodes,
region.astype('int64'))

idc = np.array(idc[:ncons, :])
c = np.array(c[:ncons, :])
B = np.zeros(c.shape[0])
self.cg = (c,idc,B)
return self.cg[0], self.cg[1], self.cg[2]

def get_elements(self):
"""
Expand Down Expand Up @@ -323,7 +269,7 @@ def get_element_gradients(self, elements = None):
# get cell corners
xi, yi, zi = self.cell_corner_indexes(c_xi, c_yi, c_zi) # global_index_to_node_index(gi)
# convert to node locations
nodes = self.node_indexes_to_position(xi, yi, zi).T
nodes = np.array(self.node_indexes_to_position(xi, yi, zi)).T

points = np.zeros((5, 4, self.n_cells, 3))
points[:, :, even_mask, :] = nodes[:, even_mask, :][self.tetra_mask_even, :, :]
Expand Down
131 changes: 90 additions & 41 deletions LoopStructural/interpolators/unstructured_tetra.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,16 +33,26 @@ def __init__(self, nodes, elements, neighbours,aabb_nsteps=None):
force nsteps for aabb, by default None
"""
self.nodes = np.array(nodes)
self.n_nodes = self.nodes.shape[0]
self.neighbours = np.array(neighbours,dtype=np.int64)
self.elements = np.array(elements,dtype=np.int64)
self.barycentre = np.sum(self.nodes[self.elements][:, :, :],
axis=1) / 4.
self.minimum = np.min(self.nodes,axis=0)
self.maximum = np.max(self.nodes,axis=0)
length = (self.maximum-self.minimum)
self.minimum -= length
self.maximum += length

self.minimum -= length*.1
self.maximum += length*.1
if aabb_nsteps==None:
box_vol = np.product(self.maximum-self.minimum)
element_volume = box_vol / (len(self.elements)/20)
# calculate the step vector of a regular cube
step_vector = np.zeros(3)
step_vector[:] = element_volume ** (1. / 3.)
# number of steps is the length of the box / step vector
aabb_nsteps = np.ceil((self.maximum-self.minimum) / step_vector).astype(int)
#make sure there is at least one cell in every dimension
aabb_nsteps[aabb_nsteps<2] =2
aabb_nsteps = np.array(aabb_nsteps,dtype=int)
step_vector = (self.maximum-self.minimum)/(aabb_nsteps-1)
self.aabb_grid = BaseStructuredSupport(self.minimum,nsteps=aabb_nsteps,step_vector=step_vector)
Expand All @@ -51,45 +61,75 @@ def __init__(self, nodes, elements, neighbours,aabb_nsteps=None):
# at the expense of speed
self.aabb_table = np.zeros((self.aabb_grid.n_elements,len(self.elements)),dtype=bool)
self.aabb_table[:] = False
self.initialise_aabb()
self._initialise_aabb()

def _initialise_aabb(self):
"""assigns the tetras to the grid cells where the bounding box
of the tetra element overlaps the grid cell. Note this is an approximation
and also does not work when the tetra bounding box is bigger than the aabb
grid cell.
of the tetra element overlaps the grid cell.
It could be changed to use the separating axis theorem, however this would require
significantly more calculations... #TODO test timing
significantly more calculations. (12 more I think).. #TODO test timing
"""
# calculate the bounding box for all tetraherdon in the mesh
# find the min/max extents for xyz
tetra_bb = np.zeros((self.elements.shape[0],8,3))
tetra_bb = np.zeros((self.elements.shape[0],19,3))
minx = np.min(self.nodes[self.elements,0],axis=1)
maxx = np.max(self.nodes[self.elements,0],axis=1)
miny = np.min(self.nodes[self.elements,1],axis=1)
maxy = np.max(self.nodes[self.elements,1],axis=1)
minz = np.min(self.nodes[self.elements,2],axis=1)
maxz = np.max(self.nodes[self.elements,2],axis=1)

# add the corners of the cube
tetra_bb[:,0,:] = np.array([minx,miny,minz]).T
tetra_bb[:,1,:] = np.array([maxx,miny,minz]).T
tetra_bb[:,2,:] = np.array([maxx,maxy,minz]).T
tetra_bb[:,3,:] = np.array([minx,maxy,minz]).T
tetra_bb[:,4,:] = np.array([minx,miny,maxz]).T
tetra_bb[:,5,:] = np.array([maxx,miny,minz]).T
tetra_bb[:,6,:] = np.array([maxx,maxy,maxz]).T
tetra_bb[:,7,:] = np.array([minx,maxy,maxz]).T

# find which
cell_index_ijk = np.array(self.aabb_grid.position_to_cell_index(tetra_bb.reshape((-1,3)))).swapaxes(0,1)
cell_index_global = cell_index_ijk[:, 0] + self.aabb_grid.nsteps_cells[ None, 0] \
* cell_index_ijk[:, 1] + self.aabb_grid.nsteps_cells[ None, 0] * \
self.aabb_grid.nsteps_cells[ None, 1] * cell_index_ijk[:, 2]
bbcorners_grid_cell = cell_index_global.reshape((tetra_bb.shape[0],tetra_bb.shape[1]))
tetra_index = np.arange(self.elements.shape[0],dtype=int)
tetra_index = np.tile(tetra_index,(8,1)).T
self.aabb_table[bbcorners_grid_cell,tetra_index] = True
ix,iy,iz = self.aabb_grid.global_index_to_cell_index(np.arange(self.aabb_grid.n_elements))
cix,ciy,ciz = self.aabb_grid.cell_corner_indexes(ix,iy,iz)
px,py,pz = self.aabb_grid.node_indexes_to_position(cix,ciy,ciz)
x_boundary = px[:,[0,1]]
y_boundary = py[:,[0,2]]
z_boundary = pz[:,[0,3]]
a = np.logical_and(minx[None,:] > x_boundary[:,None,0],minx[None,:] < x_boundary[:,None,1]) # min point between cell
b = np.logical_and(maxx[None,:] < x_boundary[:,None,1],maxx[None,:] > x_boundary[:,None,0]) # max point between cell
c = np.logical_and(minx[None,:] < x_boundary[:,None,0],maxx[None,:] > x_boundary[:,None,0]) # min point < than cell & max point > cell

x_logic = np.logical_or(np.logical_or(a,b),c)

a = np.logical_and(miny[None,:] > y_boundary[:,None,0],miny[None,:] < y_boundary[:,None,1]) # min point between cell
b = np.logical_and(maxy[None,:] < y_boundary[:,None,1],maxy[None,:] > y_boundary[:,None,0]) # max point between cell
c = np.logical_and(miny[None,:] < y_boundary[:,None,0],maxy[None,:] > y_boundary[:,None,0]) # min point < than cell & max point > cell

y_logic = np.logical_or(np.logical_or(a,b),c)


a = np.logical_and(minz[None,:] > z_boundary[:,None,0],minz[None,:] < z_boundary[:,None,1]) # min point between cell
b = np.logical_and(maxz[None,:] < z_boundary[:,None,1],maxz[None,:] > z_boundary[:,None,0]) # max point between cell
c = np.logical_and(minz[None,:] < z_boundary[:,None,0],maxz[None,:] > z_boundary[:,None,0]) # min point < than cell & max point > cell

z_logic = np.logical_or(np.logical_or(a,b),c)
logic = np.logical_and(x_logic,y_logic)
logic = np.logical_and(logic,z_logic)
# inside_x =
# # add the corners of the cube
# tetra_bb[:,0,:] = np.array([minx,miny,minz]).T
# tetra_bb[:,1,:] = np.array([maxx,miny,minz]).T
# tetra_bb[:,2,:] = np.array([maxx,maxy,minz]).T
# tetra_bb[:,3,:] = np.array([minx,maxy,minz]).T
# tetra_bb[:,4,:] = np.array([minx,miny,maxz]).T
# tetra_bb[:,5,:] = np.array([maxx,miny,minz]).T
# tetra_bb[:,6,:] = np.array([maxx,maxy,maxz]).T
# tetra_bb[:,7,:] = np.array([minx,maxy,maxz]).T
# # add centre
# tetra_bb[:,8,:] = np.array([(maxx-minx)/2,(maxy-miny)/2,(maxy-miny)/2]).T




# # find which
# cell_index_ijk = np.array(self.aabb_grid.position_to_cell_index(tetra_bb.reshape((-1,3)))).swapaxes(0,1)
# cell_index_global = cell_index_ijk[:, 0] + self.aabb_grid.nsteps_cells[ None, 0] \
# * cell_index_ijk[:, 1] + self.aabb_grid.nsteps_cells[ None, 0] * \
# self.aabb_grid.nsteps_cells[ None, 1] * cell_index_ijk[:, 2]
# bbcorners_grid_cell = cell_index_global.reshape((tetra_bb.shape[0],tetra_bb.shape[1]))
# tetra_index = np.arange(self.elements.shape[0],dtype=int)
# tetra_index = np.tile(tetra_index,(9,1)).T
self.aabb_table = logic
#[bbcorners_grid_cell,tetra_index] = True

@property
def ntetra(self):
Expand Down Expand Up @@ -175,6 +215,8 @@ def inside(self, pos):
self.step_vector[None, i] * self.nsteps_cells[None, i]
return inside

def get_elements(self):
return self.elements
def get_tetra_for_location(self, points):
"""
Determine the tetrahedron from a numpy array of points
Expand Down Expand Up @@ -219,9 +261,16 @@ def get_tetra_for_location(self, points):
c[:, 3] = vd / v
# inside = np.ones(c.shape[0],dtype=bool)
mask = np.all(c>=0,axis=1)

return vertices[mask,:,:], c[mask,:], col[mask], np.ones(col[mask].shape,dtype=bool)

verts = np.zeros((points.shape[0],4,3))
bc = np.zeros((points.shape[0],4))
tetras = np.zeros(points.shape[0],dtype='int64')
inside = np.zeros(points.shape[0],dtype=bool)

verts[row[mask],:,:] = vertices[mask,:,:]
bc[row[mask],:] = c[mask,:]
tetras[row[mask]] = col[mask]
inside[row[mask]]=True
return verts, bc, self.elements[tetras], inside


def get_element_gradients(self, elements = None):
Expand All @@ -236,16 +285,16 @@ def get_element_gradients(self, elements = None):
-------
"""
points = np.zeros((5, 4, self.n_cells, 3))
points[:, :, even_mask, :] = nodes[:, even_mask, :][self.tetra_mask_even, :, :]
points[:, :, ~even_mask, :] = nodes[:, ~even_mask, :][self.tetra_mask, :, :]

# changing order to points, tetra, nodes, coord
points = points.swapaxes(0, 2)
points = points.swapaxes(1, 2)
# points = np.zeros((5, 4, self.n_cells, 3))
# points[:, :, even_mask, :] = nodes[:, even_mask, :][self.tetra_mask_even, :, :]
# points[:, :, ~even_mask, :] = nodes[:, ~even_mask, :][self.tetra_mask, :, :]

ps = points.reshape(points.shape[0] * points.shape[1], points.shape[2], points.shape[3])
# # changing order to points, tetra, nodes, coord
# points = points.swapaxes(0, 2)
# points = points.swapaxes(1, 2)

ps = self.nodes[self.elements,:]#points.reshape(points.shape[0] * points.shape[1], points.shape[2], points.shape[3])
#vertices = self.nodes[self.elements[col,:]]
m = np.array(
[[(ps[:, 1, 0] - ps[:, 0, 0]), (ps[:, 1, 1] - ps[:, 0, 1]),
(ps[:, 1, 2] - ps[:, 0, 2])],
Expand Down Expand Up @@ -307,7 +356,7 @@ def get_neighbours(self):
-------
"""
self.neighbours
return self.neighbours



0 comments on commit f5a5f9b

Please sign in to comment.