Source code for plot3d.block

import numpy as np 
import math 
from tqdm import trange
from typing import List

[docs] class Block: """Plot3D Block definition """ def __init__(self, X:np.ndarray,Y:np.ndarray,Z:np.ndarray): """Initializes the block using all the X,Y,Z coordinates of the block Args: X (np.ndarray): All the X coordinates (i,j,k) Y (np.ndarray): All the Y coordinates (i,j,k) Z (np.ndarray): All the Z coordinates (i,j,k) """ self.IMAX,self.JMAX,self.KMAX = X.shape; self.X = X self.Y = Y self.Z = Z # Centroid self.cx = np.mean(X) self.cy = np.mean(Y) self.cz = np.mean(Z)
[docs] def scale(self,factor:float): """Scales a mesh by a certain factor Args: factor (float): _description_ """ self.X *= factor self.Y *= factor self.Z *= factor
[docs] def shift(self,shift_amount:float,direction:str="z"): """shifts the blocks by a certain amount Args: shift_amount (float): _description_ direction (str, optional): _description_. Defaults to "z". """ if direction.lower() == 'z': self.Z +=shift_amount elif direction.lower() == 'y': self.Y +=shift_amount elif direction.lower() == 'x': self.X +=shift_amount
[docs] def cylindrical(self): """Converts the block to cylindrical coordinate system. The rotation axis is assumed to be "x" direction """ self.r = np.sqrt(self.Z*self.Z + self.Y*self.Y) self.theta = np.arctan2(self.Y,self.Z)
[docs] def cell_volumes(self): """Compute volume of all cells Returns: numpy.ndarray: volume of all cells Reference: Davies, D.E. and Salmond, D.J., "Calculation of the Volume of a General Hexahedron for Flow Predicitons", AIAA Journal, vol. 23, No. 6, pp. 954-956, June 1985. It is (supposedly) exact for a hexahedral whose faces are bi-linear surfaces (i.e., the simplest surface that can be fit through the four nodes defining the face). """ X = self.X Y = self.Y Z = self.Z a = [np.zeros(shape=(self.IMAX,self.JMAX,self.KMAX)) for _ in range(9)] # face csi=const for k in range(1,self.KMAX): for j in range(1,self.JMAX): for i in range(self.IMAX): # csi-const dx1 = X[i,j,k-1] - X[i,j-1,k] dy1 = Y[i,j,k-1] - Y[i,j-1,k] dz1 = Z[i,j,k-1] - Z[i,j-1,k] dx2 = X[i,j,k] - X[i,j-1,k-1] dy2 = Y[i,j,k] - Y[i,j-1,k-1] dz2 = Z[i,j,k] - Z[i,j-1,k-1] ax = dy1*dz2 - dz1*dy2 ay = dz1*dx2 - dx1*dz2 az = dx1*dy2 - dy1*dx2 a[0][i,j,k] = ax*0.5 a[1][i,j,k] = ay*0.5 a[2][i,j,k] = az*0.5 for k in range(1,self.KMAX): for j in range(self.JMAX): # face eta=const for i in range(1,self.IMAX): dx1 = X[i,j,k] - X[i-1,j,k-1] dy1 = Y[i,j,k] - Y[i-1,j,k-1] dz1 = Z[i,j,k] - Z[i-1,j,k-1] dx2 = X[i,j,k-1] - X[i-1,j,k] dy2 = Y[i,j,k-1] - Y[i-1,j,k] dz2 = Z[i,j,k-1] - Z[i-1,j,k] ax = dy1*dz2 - dz1*dy2 ay = dz1*dx2 - dx1*dz2 az = dx1*dy2 - dy1*dx2 a[3][i,j,k] = ax*0.5 a[4][i,j,k] = ay*0.5 a[5][i,j,k] = az*0.5 # face zit=const for k in range(self.KMAX): # zit=const for j in range(1,self.JMAX): for i in range(1,self.IMAX): dx1 = X[i,j,k] - X[i-1,j-1,k] dy1 = Y[i,j,k] - Y[i-1,j-1,k] dz1 = Z[i,j,k] - Z[i-1,j-1,k] dx2 = X[i-1,j,k] - X[i,j-1,k] dy2 = Y[i-1,j,k] - Y[i,j-1,k] dz2 = Z[i-1,j,k] - Z[i,j-1,k] ax = dy1*dz2 - dz1*dy2 ay = dz1*dx2 - dx1*dz2 az = dx1*dy2 - dy1*dx2 a[6][i,j,k] = ax*0.5 a[7][i,j,k] = ay*0.5 a[8][i,j,k] = az*0.5 cf = np.zeros(shape=(6,3)) v = np.zeros(shape=(self.IMAX,self.JMAX,self.KMAX)) for k in trange(1,self.KMAX,desc='Calculating the volumes'): for j in range(1,self.JMAX): for i in range(1,self.IMAX): cf[0,0] = X[i-1,j-1,k-1] + X[i-1,j-1,k] + X[i-1,j,k-1] + X[i-1,j,k] cf[0,1] = Y[i-1,j-1,k-1] + Y[i-1,j-1,k] + Y[i-1,j,k-1] + Y[i-1,j,k] cf[0,2] = Z[i-1,j-1,k-1] + Z[i-1,j-1,k] + Z[i-1,j,k-1] + Z[i-1,j,k] cf[1,0] = X[i,j-1,k-1] + X[i,j-1,k] + X[i,j,k-1] + X[i,j,k] cf[1,1] = Y[i,j-1,k-1] + Y[i,j-1,k] + Y[i,j,k-1] + Y[i,j,k] cf[1,2] = Z[i,j-1,k-1] + Z[i,j-1,k] + Z[i,j,k-1] + Z[i,j,k] cf[2,0] = X[i-1,j-1,k-1] + X[i-1,j-1,k] + X[i,j-1,k-1] + X[i,j-1,k] cf[2,1] = Y[i-1,j-1,k-1] + Y[i-1,j-1,k] + Y[i,j-1,k-1] + Y[i,j-1,k] cf[2,2] = Z[i-1,j-1,k-1] + Z[i-1,j-1,k] + Z[i,j-1,k-1] + Z[i,j-1,k] cf[3,0] = X[i-1,j,k-1] + X[i-1,j,k] + X[i,j,k-1] + X[i,j,k] cf[3,1] = Y[i-1,j,k-1] + Y[i-1,j,k] + Y[i,j,k-1] + Y[i,j,k] cf[3,2] = Z[i-1,j,k-1] + Z[i-1,j,k] + Z[i,j,k-1] + Z[i,j,k] cf[4,0] = X[i-1,j-1,k-1] + X[i-1,j,k-1] + X[i,j-1,k-1] + X[i,j,k-1] cf[4,1] = Y[i-1,j-1,k-1] + Y[i-1,j,k-1] + Y[i,j-1,k-1] + Y[i,j,k-1] cf[4,2] = Z[i-1,j-1,k-1] + Z[i-1,j,k-1] + Z[i,j-1,k-1] + Z[i,j,k-1] cf[5,0] = X[i-1,j-1,k] + X[i-1,j,k] + X[i,j-1,k] + X[i,j,k] cf[5,1] = Y[i-1,j-1,k] + Y[i-1,j,k] + Y[i,j-1,k] + Y[i,j,k] cf[5,2] = Z[i-1,j-1,k] + Z[i-1,j,k] + Z[i,j-1,k] + Z[i,j,k] vol12=0 for n in range(0,2): # n = 0,1 for l in range(0,3): # l = 0,1,2 vol12 += math.pow(-1,n+1) *( +cf[n,l]*a[l][i-1+n,j,k] +cf[2+n,l]*a[3+l][i,j-1+n,k] +cf[4+n,l]*a[6+l][i,j,k-1+n]) v[i,j,k]= vol12/12 return v
def checkCollinearity(v1:np.ndarray, v2:np.ndarray): # Calculate their cross product cross_P = np.cross(v1,v2) # Check if their cross product # is a NULL Vector or not if (cross_P[0] == 0 and cross_P[1] == 0 and cross_P[2] == 0): return True else: return False def calculate_outward_normals(block:Block): # Calculate Normals X = block.X Y = block.Y Z = block.Z imax = block.IMAX jmax = block.JMAX kmax = block.KMAX # IMAX - Normal should be out of the page # Normals I direction: IMIN https://www.khronos.org/opengl/wiki/Calculating_a_Surface_Normal x = [X[0,0,0],X[0,jmax,0],X[0,0,kmax]] y = [Y[0,0,0],Y[0,jmax,0],Y[0,0,kmax]] z = [Z[0,0,0],Z[0,jmax,0],Z[0,0,kmax]] u = np.array([x[1]-x[0],y[1]-y[0],z[1]-z[0]]) v = np.array([x[2]-x[0],y[2]-y[0],z[2]-z[0]]) n_imin = np.cross(v1,v2) # Normals I direction: IMAX x = [X[imax,0,0],X[imax,jmax,0],X[imax,0,kmax]] y = [Y[imax,0,0],Y[imax,jmax,0],Y[imax,0,kmax]] z = [Z[imax,0,0],Z[imax,jmax,0],Z[imax,0,kmax]] v1 = np.array([x[1]-x[0],y[1]-y[0],z[1]-z[0]]) v2 = np.array([x[2]-x[0],y[2]-y[0],z[2]-z[0]]) n_imax = np.cross(v1,v2) # Normals J direction: JMIN x = [X[0,0,0],X[imax,0,0],X[0,0,kmax]] y = [Y[0,0,0],Y[imax,0,0],Y[0,0,kmax]] z = [Z[0,0,0],Z[imax,0,0],Z[0,0,kmax]] v1 = np.array([x[1]-x[0],y[1]-y[0],z[1]-z[0]]) v2 = np.array([x[2]-x[0],y[2]-y[0],z[2]-z[0]]) n_jmin = np.cross(v1,v2) # Normals J direction: JMAX x = [X[0,jmax,0],X[imax,jmax,0],X[0,jmax,kmax]] y = [Y[0,jmax,0],Y[imax,jmax,0],Y[0,jmax,kmax]] z = [Z[0,jmax,0],Z[imax,jmax,0],Z[0,jmax,kmax]] v1 = np.array([x[1]-x[0],y[1]-y[0],z[1]-z[0]]) v2 = np.array([x[2]-x[0],y[2]-y[0],z[2]-z[0]]) n_jmax = np.cross(v1,v2) # Normals K direction: KMIN x = [X[imax,0,0],X[0,jmax,0],X[0,0,0]] y = [Y[imax,0,0],Y[0,jmax,0],Y[0,0,0]] z = [Z[imax,0,0],Z[0,jmax,0],Z[0,0,0]] v1 = np.array([x[1]-x[0],y[1]-y[0],z[1]-z[0]]) v2 = np.array([x[2]-x[0],y[2]-y[0],z[2]-z[0]]) n_kmin = np.cross(v1,v2) # Normals K direction: KMAX x = [X[imax,0,kmax],X[0,jmax,kmax],X[0,0,kmax]] y = [Y[imax,0,kmax],Y[0,jmax,kmax],Y[0,0,kmax]] z = [Z[imax,0,kmax],Z[0,jmax,kmax],Z[0,0,kmax]] v1 = np.array([x[1]-x[0],y[1]-y[0],z[1]-z[0]]) v2 = np.array([x[2]-x[0],y[2]-y[0],z[2]-z[0]]) n_kmax = np.cross(v1,v2) return n_imin,n_jmin,n_kmin,n_imax,n_jmax,n_kmax
[docs] def reduce_blocks(blocks:List[Block],factor:int): """reduce the blocks by a factor of (factor) Args: blocks (List[Block]): list of blocks to reduce in size factor (int, optional): Number of indicies to skip . Defaults to 2. Returns: [type]: [description] """ for i in range(len(blocks)): blocks[i].X = blocks[i].X[::factor,::factor,::factor] blocks[i].Y = blocks[i].Y[::factor,::factor,::factor] blocks[i].Z = blocks[i].Z[::factor,::factor,::factor] blocks[i].IMAX,blocks[i].JMAX,blocks[i].KMAX = blocks[i].X.shape return blocks