Source code for plot3d.facefunctions

from typing import Dict, List
from .listfunctions import unique_pairs
from .block import Block, reduce_blocks
from .face import Face 
from copy import deepcopy
import math
import numpy as np


[docs] def get_outer_faces(block1:Block): """Get the outer faces of a block Args: block1 (Block): A plot3D block Returns: List[Face]: Non matching faces of the block List[(Face,Face)]: Matching faces inside the block """ I = [0,block1.IMAX-1] # Python index starts at 0, need to subtract 1 for it to get the i,j,k J = [0,block1.JMAX-1] K = [0,block1.KMAX-1] # Create the outer faces faces = list() face = Face(4) i=I[0] for j in J: for k in K: face.add_vertex(block1.X[i,j,k], block1.Y[i,j,k], block1.Z[i,j,k],i,j,k) faces.append(face) face = Face(4) i=I[1] for j in J: for k in K: face.add_vertex(block1.X[i,j,k], block1.Y[i,j,k], block1.Z[i,j,k],i,j,k) faces.append(face) face = Face(4) j=J[0] for i in I: for k in K: face.add_vertex(block1.X[i,j,k], block1.Y[i,j,k], block1.Z[i,j,k],i,j,k) faces.append(face) face = Face(4) j=J[1] for i in I: for k in K: face.add_vertex(block1.X[i,j,k], block1.Y[i,j,k], block1.Z[i,j,k],i,j,k) faces.append(face) face = Face(4) k=K[0] for i in I: for j in J: face.add_vertex(block1.X[i,j,k], block1.Y[i,j,k], block1.Z[i,j,k],i,j,k) faces.append(face) face = Face(4) k=K[1] for i in I: for j in J: face.add_vertex(block1.X[i,j,k], block1.Y[i,j,k], block1.Z[i,j,k],i,j,k) faces.append(face) # Check if faces match each other matching = list() non_matching = list() for i in range(len(faces)): matchFound = False for j in range(len(faces)): if (i!=j and faces[i].vertices_equals(faces[j])): matching.append((i,j)) matchFound = True if not matchFound: non_matching.append(faces[i]) # these are guaranteed to be exterior matching = list(unique_pairs(matching)) matching = [(faces[i],faces[j]) for i,j in matching] # Make sure normals do not intersect # block_center_to_face_center = block1.cx return non_matching, matching # these should be the outer faces
[docs] def create_face_from_diagonals(block:Block,imin:int,jmin:int,kmin:int,imax:int,jmax:int,kmax:int): """Creates a face on a block given a the diagonals defined as (IMIN,JMIN,KMIN), (IMAX, JMAX, KMAX) Args: block (Block): Block to create a face on imin (int): Lower Corner IMIN jmin (int): Lower Corner JMIN kmin (int): Lower Corner KMIN imax (int): Upper Corner IMAX jmax (int): Upper Corner JMAX kmax (int): Upper Corner Returns: (Face): Face created from diagonals """ newFace = Face(4) # This is because two of the corners either imin or imax can be equal if imin==imax: i = imin for j in [jmin,jmax]: for k in [kmin,kmax]: x = block.X[i,j,k] y = block.Y[i,j,k] z = block.Z[i,j,k] newFace.add_vertex(x,y,z,i,j,k) elif jmin==jmax: j = jmin for i in [imin,imax]: for k in [kmin,kmax]: x = block.X[i,j,k] y = block.Y[i,j,k] z = block.Z[i,j,k] newFace.add_vertex(x,y,z,i,j,k) elif kmin==kmax: k = kmin for i in [imin,imax]: for j in [jmin,jmax]: x = block.X[i,j,k] y = block.Y[i,j,k] z = block.Z[i,j,k] newFace.add_vertex(x,y,z,i,j,k) return newFace
[docs] def find_connected_faces(face_to_search:Face,outer_faces:List[Face],connectivity_matrix:np.ndarray,blocks:List[Block]): """Recursive program to return all the matching faces. Note faces must have the same I,J,K definition so faces will be matching if for example: Face1 IMIN=IMAX and Face2 IMIN=IMAX and they share a common edge (2 vertices) Args: face_to_search (Face): This is the face to search for outer_faces (List[Face]): List of outer faces connectivity_matrix (np.ndarray): block connectivity matrix searched_faces (List[Face]): Leave this as blank Returns: List[Face]: list of all faces that connect with face_to_search and it's neighbors. Beware of duplicates. """ connectivity_matrix = connectivity_matrix.copy() all_matching_faces = [face_to_search] faces_to_search = [face_to_search] faces_searched = [] while len(faces_to_search)>0: matching_faces = list() for face_to_search in faces_to_search: n1 = face_to_search.normal(blocks[face_to_search.BlockIndex]) selected_block_indx = face_to_search.BlockIndex connected_block_indices = np.where(connectivity_matrix[selected_block_indx,:]==1)[0] faces_to_check = [o for o in outer_faces if o.BlockIndex in connected_block_indices.tolist()] for f in faces_to_check: if (len(face_to_search.match_indices(f))==2): n2 = f.normal(blocks[f.BlockIndex]) angle = abs(math.degrees(math.acos(np.dot(n1,n2)/(np.linalg.norm(n1)*np.linalg.norm(n2))))) if angle>90: angle = 180-angle if angle<30: connectivity_matrix[selected_block_indx, f.BlockIndex] = 0 connectivity_matrix[f.BlockIndex, selected_block_indx] = 0 matching_faces.append(f) faces_searched.extend(faces_to_search) matching_faces = [m for m in matching_faces if m not in all_matching_faces] faces_to_search.clear() faces_to_search.extend(matching_faces) all_matching_faces.extend(matching_faces) matching_faces.clear() all_matching_faces = list(set(all_matching_faces)) return all_matching_faces
[docs] def find_closest_block(blocks:List[Block],x:np.ndarray,y:np.ndarray,z:np.ndarray,centroid:np.ndarray,translational_direction:str="x",minvalue:bool=True): """Find the closest block to an extreme in the x,y, or z direction and returns the targetting point. Target point is the reference point where we want the closest block and the closest face Args: x (np.ndarray): x coordinate of all the blocks' centroid y (np.ndarray): y coordinate of all the blocks' centroid z (np.ndarray): z coordinate of all the blocks' centroid centroid (np.ndarray): centroid (cx,cy,cz) translational_direction (str, optional): _description_. Defaults to "x". minvalue (bool, optional): _description_. Defaults to True. Returns: (tuple): containing *selected block index* (int): index of closest block *target_x* (float): this is the x value where selected block is closest to *target_y* (float): this is the y value where selected block is closest to *target_z* (float): this is the z value where selected block is closest to """ cx = centroid[0]; cy = centroid[1]; cz = centroid[2] target_x = cx; target_y = cy; target_z = cz if translational_direction=="x": xmins = [b.X.min() for b in blocks] xmaxes = [b.X.max() for b in blocks] xmin = min(xmins) xmax = max(xmaxes) dx = xmax - xmin if minvalue: target_x = xmin - dx*0.5 selected_block_indx = np.argmin(np.sqrt((target_x-x)**2 + (cy-y)**2 + (cz-z)**2)) else: target_x = xmax + dx*0.5 selected_block_indx = np.argmin(np.sqrt((target_x-x)**2 + (cy-y)**2 + (cz-z)**2)) target_y= cy target_z = cz elif translational_direction=="y": ymins = [b.Y.min() for b in blocks] ymaxes = [b.Y.max() for b in blocks] ymin = min(ymins) ymax = max(ymaxes) dy = ymax-ymin if minvalue: target_y = ymin - dy*0.5 selected_block_indx = np.argmin(np.sqrt((cx-x)**2 + (target_y-y)**2 + (cz-z)**2)) else: target_y = ymax + dy*0.5 selected_block_indx = np.argmin(np.sqrt((cx-x)**2 + (target_y-y)**2 + (cz-z)**2)) target_x = cx target_z = cz else: # translational_direction=="z": zmins = [b.Z.min() for b in blocks] zmaxes = [b.Z.max() for b in blocks] zmin = min(zmins) zmax = max(zmaxes) dz = zmax - zmin if minvalue: target_z = zmin - dz*0.5 selected_block_indx = np.argmin(np.sqrt((cx-x)**2 + (cy-y)**2 + (target_z-z)**2)) else: target_z = zmax + dz*0.5 selected_block_indx = np.argmin(np.sqrt((cx-x)**2 + (cy-y)**2 + (target_z-z)**2)) target_x = cx target_y = cy return selected_block_indx,target_x,target_y,target_z
[docs] def find_bounding_faces(blocks:List[Block],connectivity_matrix:np.ndarray,outer_faces:List[Dict[str,int]], direction:str="z"): """Finds bounding faces in x,y,z direction so think of this as the xmin and xmax faces of a block. This is useful for translational periodicity where you want the bounds to check Args: blocks (List[Block]): List of blocks outer_faces (List[Dict[str,int]]): outer faces as a list of dictionaries direction (str, optional): Direction to search for. Defaults to "z". Returns: (tuple) containing: *lower_connected_faces_export* (List[Dict[str,int]]): Export ready version of lower/left connected faces *upper_connected_faces_export* (List[Dict[str,int]]): Export ready version of upper/right connected faces *lower_connected_faces* (List[face]): List of lower/left connected faces *upper_connected_faces* (List[face]): List of upper/right connected faces """ gcd_array = list() # Find the gcd of all the blocks for block_indx in range(len(blocks)): block = blocks[block_indx] gcd_array.append(math.gcd(block.IMAX-1, math.gcd(block.JMAX-1, block.KMAX-1))) gcd_to_use = min(gcd_array) # You need to use the minimum gcd otherwise 1 block may not exactly match the next block. They all have to be scaled the same way. blocks = reduce_blocks(deepcopy(blocks),gcd_to_use) xyz_array = np.array([(b.cx, b.cy, b.cz) for b in blocks]); cx = np.mean(xyz_array[:,0]) # Centroid cy = np.mean(xyz_array[:,1]) cz = np.mean(xyz_array[:,2]) x = xyz_array[:,0]; y = xyz_array[:,1]; z = xyz_array[:,2] if len(outer_faces) == 0: for i,b in enumerate(blocks): outer,_ = get_outer_faces(b) [o.set_block_index(i) for o in outer] outer_faces.extend(outer) outer_faces_all = outer_faces else: outer_faces_all = outer_face_dict_to_list(blocks,outer_faces,gcd_to_use) # Find closest face to the min value selected_block_index,tx,ty,tz = find_closest_block(blocks,x,y,z,np.array([cx,cy,cz]),direction,minvalue=True) faces = [f for f in outer_faces_all if f.BlockIndex == selected_block_index] # Pick all the faces for the selected block min_face_indx = np.argmin(np.array([np.sqrt((f.cx-tx)**2+(f.cy-ty)**2+(f.cz-tz)**2) for f in faces])) min_face = faces[min_face_indx] # Just need this selected_block_index,tx,ty,tz = find_closest_block(blocks,x,y,z,np.array([cx,cy,cz]),direction,minvalue=False) faces = [f for f in outer_faces_all if f.BlockIndex == selected_block_index] # Pick all the faces for the selected block max_face_indx = np.argmin(np.array([np.sqrt((f.cx-tx)**2+(f.cy-ty)**2+(f.cz-tz)**2) for f in faces])) max_face = faces[max_face_indx] # Also need this # Search face connectivity in block connectivity matrix connectivity_matrix = connectivity_matrix - np.eye(connectivity_matrix.shape[0],dtype=np.int8) # turn off connections with itself outer_faces_all = [o for o in outer_faces_all if o.BlockIndex != min_face.BlockIndex] outer_faces_all = [o for o in outer_faces_all if o.BlockIndex != max_face.BlockIndex] # !print("Recursively searching for connected faces") lower_connected_faces = find_connected_faces(min_face,outer_faces_all,connectivity_matrix,blocks) upper_connected_faces = find_connected_faces(max_face,outer_faces_all,connectivity_matrix,blocks) lower_connected_faces.append(min_face) upper_connected_faces.append(max_face) # Unscale faces lower_connected_faces = list(set(lower_connected_faces)) upper_connected_faces = list(set(upper_connected_faces)) for l in lower_connected_faces: l.I*=gcd_to_use l.J*=gcd_to_use l.K*=gcd_to_use for u in upper_connected_faces: u.I*=gcd_to_use u.J*=gcd_to_use u.K*=gcd_to_use lower_connected_faces_export = [l.to_dict() for l in lower_connected_faces] upper_connected_faces_export = [u.to_dict() for u in upper_connected_faces] return lower_connected_faces_export,upper_connected_faces_export, lower_connected_faces, upper_connected_faces
[docs] def split_face(face_to_split:Face, block:Block,imin:int,jmin:int,kmin:int,imax:int,jmax:int,kmax:int): """Splits a face with another face within the same block picture the split as a two rectangles inside each other Args: face_to_split (Face): Face on the block to be split block (Block): Block the split is occuring on imin (int): IMIN index of the split (diagonals) jmin (int): JMIN index of the split (diagonals) kmin (int): KMIN index of the split (diagonals) imax (int): IMAX index of the split jmax (int): JMAX index of the split kmax (int): KMAX index of the split :: left face top face right face ________ __ __ __ | __ | | | |__| | | __ | |__| | | | __ | | |__| face_to_split/center face |________| |__| |__| |__| bottom face Returns: [List[Faces]]: List of unique faces from the split """ center_face = create_face_from_diagonals(block, imin=imin,imax=imax, jmin=jmin,jmax=jmax, kmin=kmin,kmax=kmax) if kmin == kmax: # In the picture above Horizontal = i, vertical = j left_face = create_face_from_diagonals(block, imin=face_to_split.IMIN,imax=imin, jmin=face_to_split.JMIN,jmax=face_to_split.JMAX, kmin=kmin, kmax=kmax) right_face = create_face_from_diagonals(block, imin=imax, imax=face_to_split.IMAX, jmin=face_to_split.JMIN, jmax=face_to_split.JMAX, kmin=kmin, kmax=kmax) top_face = create_face_from_diagonals(block, imin=imin,imax=imax, jmin=jmax,jmax=face_to_split.JMAX, kmin=kmin,kmax=kmax) bottom_face = create_face_from_diagonals(block, imin=imin,imax=imax, jmin=face_to_split.JMIN,jmax=jmin, kmin=kmin,kmax=kmax) elif (imin==imax): # In the picture above Horizontal = j, vertical = k left_face = create_face_from_diagonals(block, imin=imin,imax=imax, jmin=face_to_split.JMIN, jmax=jmin, kmin=face_to_split.KMIN,kmax=face_to_split.KMAX) right_face = create_face_from_diagonals(block, imin=imin,imax=imax, jmin=jmax, jmax=face_to_split.JMAX, kmin=face_to_split.KMIN,kmax=face_to_split.KMAX) top_face = create_face_from_diagonals(block, imin=imin,imax=imax, jmin=jmin,jmax=jmax, kmin=kmax,kmax=face_to_split.KMAX) bottom_face = create_face_from_diagonals(block, imin=imin,imax=imax, jmin=jmin,jmax=jmax, kmin=face_to_split.KMIN,kmax=kmin) elif (jmin==jmax): # In the picture above Horizontal = i, vertical = k left_face = create_face_from_diagonals(block, imin=face_to_split.IMIN,imax=imin, jmin=jmin,jmax=jmax, kmin=face_to_split.KMIN,kmax=face_to_split.KMAX) right_face = create_face_from_diagonals(block, imin=imax,imax=face_to_split.IMAX, jmin=jmin,jmax=jmax, kmin=face_to_split.KMIN,kmax=face_to_split.KMAX) top_face = create_face_from_diagonals(block, imin=imin,imax=imax, jmin=jmin,jmax=jmax, kmin=kmax,kmax=face_to_split.KMAX) bottom_face = create_face_from_diagonals(block, imin=imin, imax=imax, jmin=jmin, jmax=jmax, kmin=face_to_split.KMIN, kmax=kmin) faces = [top_face,bottom_face,left_face,right_face] faces = [f for f in faces if not f.isEdge and not f.index_equals(center_face)] # Remove edges [f.set_block_index(face_to_split.blockIndex) for f in faces] return faces
[docs] def find_face_nearest_point(faces:List[Face], x:float,y:float,z:float): """Find a face nearest to a given point Args: blocks (List[Block]): List of blocks faces (List[Face]): List of faces x (float): x coordinate of a reference point y (float): y coordinate of a reference point z (float): z coordinate of a reference point """ n = list(range(len(faces))) dv = list() for i in n: dx = x-faces[i].cx dy = y-faces[i].cy dz = z-faces[i].cz dv.append(math.sqrt(dx*dx + dy*dy + dz*dz)) face_index = np.argmin(np.array(dv)) return face_index
[docs] def outer_face_dict_to_list(blocks:List[Block],outer_faces:List[Dict[str,int]],gcd:int=1) -> List[Face]: """Converts a list of dictionary face representations to a list of faces. Use this only for outer faces Args: blocks (List[Block]): List of blocks outer_faces (List[Dict[str,int]]): List of outer faces represented as a dictionary gcd (int, optional): Greatst common divisor. Defaults to 1. Returns: List[Face]: List of Face objects """ outer_faces_all = list() for o in outer_faces: face = create_face_from_diagonals(blocks[o['block_index']], int(o['IMIN']/gcd), int(o['JMIN']/gcd), int(o['KMIN']/gcd), int(o['IMAX']/gcd), int(o['JMAX']/gcd), int(o['KMAX']/gcd)) if 'id' in o.keys(): face.id = o['id'] face.set_block_index(o['block_index']) outer_faces_all.append(face) return outer_faces_all
[docs] def match_faces_dict_to_list(blocks:List[Block],matched_faces:List[Dict[str,int]],gcd:int=1): """Converts a list of dictionaries representing matched faces to a list of Faces Args: blocks (List[Block]): List of blocks matched_faces (List[Dict[str,int]]): List of matched faces represented as a dictionary gcd (int, optional): _description_. Defaults to 1. Returns: _type_: _description_ """ matched_faces_all = list() for _,m in enumerate(matched_faces): face1 = create_face_from_diagonals(blocks[m['block1']['block_index']], int(m['block1']['IMIN']/gcd), int(m['block1']['JMIN']/gcd), int(m['block1']['KMIN']/gcd), int(m['block1']['IMAX']/gcd), int(m['block1']['JMAX']/gcd), int(m['block1']['KMAX']/gcd)) face2 = create_face_from_diagonals(blocks[m['block2']['block_index']], int(m['block2']['IMIN']/gcd), int(m['block2']['JMIN']/gcd), int(m['block2']['KMIN']/gcd), int(m['block2']['IMAX']/gcd), int(m['block2']['JMAX']/gcd), int(m['block2']['KMAX']/gcd)) face1.set_block_index(m['block1']['block_index']) if 'id' in m['block1'].keys(): face1.id = m['block1']['id'] face2.set_block_index(m['block2']['block_index']) if 'id' in m['block2'].keys(): face2.id = m['block2']['id'] matched_faces_all.append(face1) matched_faces_all.append(face2) return matched_faces_all