diff --git a/python/examples/block_test/block_test.py b/python/examples/block_test/block_test.py index fc733e8..30d5d5a 100644 --- a/python/examples/block_test/block_test.py +++ b/python/examples/block_test/block_test.py @@ -1,25 +1,77 @@ -import os, sys +import sys, os, pickle +import numpy as np sys.path.insert(0,'../../') -sys.path.insert(0,'../Cascade') -from plot3d import write_plot3D, read_plot3D, periodicity -from plot3d import find_matching_blocks, get_outer_faces, connectivity, connectivity_fast -from glennht_con import export_to_glennht_conn -import pickle - - -if not os.path.exists('connectivity.pickle'): - blocks = read_plot3D('connectivity_test.x', binary = False) - # Block 1 is the blade O-Mesh k=0 - # outer_faces, _ = get_outer_faces(blocks[0]) # lets check - face_matches, outer_faces_formatted = connectivity_fast(blocks) - # face_matches2, outer_faces_formatted2 = connectivity(blocks) - with open('connectivity.pickle','wb') as f: - [m.pop('match',None) for m in face_matches] # Remove the dataframe - pickle.dump({"face_matches":face_matches, "outer_faces":outer_faces_formatted},f) - -with open('connectivity.pickle','rb') as f: - data = pickle.load(f) - face_matches = data['face_matches'] - outer_faces = data['outer_faces'] - -export_to_glennht_conn(face_matches,outer_faces,'block_test') +from plot3d import read_plot3D, connectivity_fast,translational_periodicity, write_plot3D, Direction, split_blocks, block_connection_matrix, find_bounding_faces +from plot3d import outer_face_dict_to_list, match_faces_dict_to_list + +#%% Find connectivity +def dump_data(data): + with open('block_data.pickle','wb') as f: + pickle.dump(data,f) + +def read_data(): + with open('block_data.pickle','rb') as f: + return pickle.load(f) + +blocks = read_plot3D('iso65_64blocks.xyz',True) + +if not os.path.exists(f'block_data.pickle'): + print('Finding connectivity') + face_matches, outer_faces = connectivity_fast(blocks) + [m.pop('match',None) for m in face_matches] # Remove the dataframe + print('Organizing split and outerfaces') + all_faces = match_faces_dict_to_list(blocks,face_matches) + all_faces.extend(outer_face_dict_to_list(blocks,outer_faces)) + all_faces = [m.to_dict() for m in all_faces] + data = { + "face_matches":face_matches, + "outer_faces":outer_faces, + "all_faces":all_faces + } + dump_data(data) + print('Creating block connection matrix') + c = block_connection_matrix(blocks,all_faces) + data["connectivity_matrix"]=c + dump_data(data) + +data = read_data() +all_faces = data['all_faces'] +connectivity_matrix = data['connectivity_matrix'] + +#%% Find bounding Faces +lower_bound, upper_bound,_,_ = find_bounding_faces(blocks,connectivity_matrix,all_faces,"z") +left_bound, right_bound,_,_ = find_bounding_faces(blocks,connectivity_matrix,all_faces,"y") +data['lower_bound'] = lower_bound +data['upper_bound'] = upper_bound +data['left_bound'] = left_bound +data['right_bound'] = right_bound +dump_data(data) + +#%% Use bounding faces to find periodicity +data = read_data() +lower_bound = data['lower_bound']; upper_bound = data['upper_bound'] +left_bound = data['left_bound']; right_bound = data['right_bound'] +x_periodic_faces_export, periodic_faces = translational_periodicity(blocks,left_bound,right_bound,translational_direction='x') +y_periodic_faces_export, periodic_faces = translational_periodicity(blocks,left_bound,right_bound,translational_direction='y') +z_periodic_faces_export, periodic_faces = translational_periodicity(blocks,lower_bound,upper_bound,translational_direction='z') +data['x_periodic'] = z_periodic_faces_export +data['z_periodic'] = z_periodic_faces_export +data['y_periodic'] = y_periodic_faces_export +dump_data(data) + +#%% Lets check with faces are not periodic in the y-direction +data = read_data() +y_periodic_faces_export = data['y_periodic'] +left_bound = data['left_bound']; right_bound = data['right_bound'] + +left_periodic_blocks_found = [p['block1']['block_index'] for p in y_periodic_faces_export] +left_faces_missing = [l for l in left_bound if l['block_index'] not in left_periodic_blocks_found] +print('Left faces missing') +[print(l) for l in left_faces_missing] + +right_periodic_blocks_found = [p['block2']['block_index'] for p in y_periodic_faces_export] +right_faces_missing = [r for r in right_bound if r['block_index'] not in right_periodic_blocks_found] +print('Right faces missing') +[print(r) for r in right_faces_missing] + +# Lets find out why it's missing \ No newline at end of file diff --git a/python/examples/block_test/block_tutorial.ipynb b/python/examples/block_test/block_tutorial.ipynb new file mode 100644 index 0000000..8ba0f82 --- /dev/null +++ b/python/examples/block_test/block_tutorial.ipynb @@ -0,0 +1,393 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Translated Periodicity using a simple Block\n", + "The mesh shown below has periodic faces on the z and y axis. The code below shows how we find the translated periodicity in the y and z directions.\n", + "\n", + "Steps Overview\n", + "1. Find Connectivity - This splits the faces of all the blocks so that they match \n", + "2. Build the Block connectivity matrix using connectivity - Uses the split faces to determine which blocks are connected. value of 1 = connected, 0 = not searched yet, -1 = not connected \n", + "3. Search for connected faces on the z and y axis - Uses the block connectivity to find the outer connected faces. So, all the top faces for instance. \n", + "4. Find translated periodicity - Uses the connected faces e.g. the left and right, or top and bottom, or front and back. Checks the connected faces to see if any of them match faces from the other side. \n", + "\n", + "![image.png](tutorial_images/cmc009-overview.png)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Install Plot3D\n", + "!pip install plot3d" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "# Import functions\n", + "import os, pickle\n", + "from plot3d import read_plot3D, connectivity_fast,translational_periodicity, write_plot3D, Direction, split_blocks, block_connection_matrix, find_bounding_faces\n", + "from plot3d import outer_face_dict_to_list, match_faces_dict_to_list\n", + "import numpy as np\n", + "#%% Find connectivity \n", + "def dump_data(data):\n", + " with open('block_data.pickle','wb') as f:\n", + " pickle.dump(data,f)\n", + "\n", + "def read_data():\n", + " with open('block_data.pickle','rb') as f:\n", + " return pickle.load(f)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of blocks: 64\n" + ] + } + ], + "source": [ + "# Load the mesh\n", + "blocks = read_plot3D('iso65_64blocks.p3d',False)\n", + "write_plot3D('iso65_64blocks.xyz',blocks,True)\n", + "n = len(blocks)\n", + "print(f'Number of blocks: {n}')" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 1: Find Connectivity\n", + "The code below finds the connectivity and creates the pickle file storing the connectivity " + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Finding connectivity\n", + "gcd to use 16\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Finding nearest blocks comparing 62 with 63: 100%|██████████| 2016/2016 [00:05<00:00, 382.76it/s]\n", + "Checking connections block 62 with 63: 100%|██████████| 363/363 [00:18<00:00, 19.73it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Organizing split and outerfaces\n" + ] + } + ], + "source": [ + "print('Finding connectivity')\n", + "face_matches, outer_faces = connectivity_fast(blocks)\n", + "[m.pop('match',None) for m in face_matches] # Remove the dataframe to save as pickle file \n", + "print('Organizing split and outerfaces')\n", + "all_faces = match_faces_dict_to_list(blocks,face_matches)\n", + "all_faces.extend(outer_face_dict_to_list(blocks,outer_faces))\n", + "all_faces = [m.to_dict() for m in all_faces]\n", + "data = {\n", + " \"face_matches\":face_matches, \n", + " \"outer_faces\":outer_faces,\n", + " \"all_faces\":all_faces\n", + " }\n", + "dump_data(data)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 2: Build Block to Block connectivity matrix\n", + "Block to block connectivity makes it easier to search for which blocks are connected. For example `print(c[10,4])` if this is 1 then that means blocks 10 and 4 are touching. \n", + "\n", + "Figure below shows an example of the connections for block 704. \n", + "\n", + "![block connectivity](tutorial_images/cmc009-block-connectivity.png)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Creating block connection matrix\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Building block to block connectivity matrix: checking 62: 100%|██████████| 2016/2016 [00:02<00:00, 999.50it/s] \n" + ] + } + ], + "source": [ + "data = read_data()\n", + "print('Creating block connection matrix')\n", + "c = block_connection_matrix(blocks,all_faces)\n", + "data[\"connectivity_matrix\"]=c\n", + "dump_data(data)\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These are the blocks where block 2 is connected to. We have blocks 0,3,6, and 16. Block 2 is automtically connected with itself. " + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "connectivity matrix shape: 64 by 64\n" + ] + }, + { + "data": { + "text/plain": [ + "(array([ 0, 2, 3, 6, 16], dtype=int64),)" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "print(\"connectivity matrix shape: \" + str(c.shape[0]) + ' by ' + str(c.shape[1]))\n", + "np.where(c[2,:]==1)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 3: Find bounding faces\n", + "These are faces on either the left or right side of the block." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "data = read_data() \n", + "all_faces = data['all_faces']\n", + "connectivity_matrix = data['connectivity_matrix']\n", + "\n", + "#%% Find bounding Faces\n", + "lower_bound, upper_bound,_,_ = find_bounding_faces(blocks,connectivity_matrix,all_faces,\"z\")\n", + "left_bound, right_bound,_,_ = find_bounding_faces(blocks,connectivity_matrix,all_faces,\"y\")\n", + "data['lower_bound'] = lower_bound\n", + "data['upper_bound'] = upper_bound\n", + "data['left_bound'] = left_bound\n", + "data['right_bound'] = right_bound\n", + "dump_data(data)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 4: Find Periodicity\n", + "In this step we determine the connected faces in the y-direction as well as the z-direction. The program will automatically inform the user if it's not periodic due to not all faces matching. If this is the case then the faces should be plotted using paraview and the user needs to check with their grid generator or manually add in the connectivity. " + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Checking connections block 17 with 54: 0%| | 0/16 [00:00 ijkmax): + temp = ijkmax + ijkmax = ijkmin + ijkmin = temp + return ijkmin, ijkmax + + for b in blocks_to_extract: # Block indicies + block_source,block_display,LUT = ExtractBlocks(plot3D_source,View,[b]) + RenameSource('Block '+str(b), block_source) + block_source = FindSource('Block '+str(b)) + + for surface_indx, o in enumerate(outer_faces): + # Add Plots for Outer Faces + if o['block_index'] == b: + voi = [o['IMIN'], o['IMAX'], o['JMIN'], o['JMAX'],o['KMIN'], o['KMAX']] + CreateSubset(block_source, voi, name='outer_face '+str(surface_indx),opacity=1) + + # Plot the periodic faces + for periodic_indx, p in enumerate(periodic_faces): + # Add Plots for Outer Faces + if p['block1']['block_index'] == b and p['block2']['block_index'] == b: # Periodicity within the block + p['block1']['IMIN'], p['block1']['IMAX'] = check_and_swap(p['block1']['IMIN'], p['block1']['IMAX']) + p['block1']['JMIN'], p['block1']['JMAX'] = check_and_swap(p['block1']['JMIN'], p['block1']['JMAX']) + p['block1']['KMIN'], p['block1']['KMAX'] = check_and_swap(p['block1']['KMIN'], p['block1']['KMAX']) + voi = [p['block1']['IMIN'], p['block1']['IMAX'], p['block1']['JMIN'], p['block1']['JMAX'],p['block1']['KMIN'], p['block1']['KMAX']] + CreateSubset(block_source, voi, name='periodic '+str(periodic_indx)) + + p['block2']['IMIN'], p['block2']['IMAX'] = check_and_swap(p['block2']['IMIN'], p['block2']['IMAX']) + p['block2']['JMIN'], p['block2']['JMAX'] = check_and_swap(p['block2']['JMIN'], p['block2']['JMAX']) + p['block2']['KMIN'], p['block2']['KMAX'] = check_and_swap(p['block2']['KMIN'], p['block2']['KMAX']) + voi = [p['block2']['IMIN'], p['block2']['IMAX'], p['block2']['JMIN'], p['block2']['JMAX'],p['block2']['KMIN'], p['block2']['KMAX']] + CreateSubset(block_source, voi, name='periodic '+str(periodic_indx)) + + elif p['block1']['block_index'] == b or p['block2']['block_index'] == b: # Periodicity from block to block + if p['block1']['block_index'] == b: + p['block1']['IMIN'], p['block1']['IMAX'] = check_and_swap(p['block1']['IMIN'], p['block1']['IMAX']) + p['block1']['JMIN'], p['block1']['JMAX'] = check_and_swap(p['block1']['JMIN'], p['block1']['JMAX']) + p['block1']['KMIN'], p['block1']['KMAX'] = check_and_swap(p['block1']['KMIN'], p['block1']['KMAX']) + voi = [p['block1']['IMIN'], p['block1']['IMAX'], p['block1']['JMIN'], p['block1']['JMAX'],p['block1']['KMIN'], p['block1']['KMAX']] + else: + p['block2']['IMIN'], p['block2']['IMAX'] = check_and_swap(p['block2']['IMIN'], p['block2']['IMAX']) + p['block2']['JMIN'], p['block2']['JMAX'] = check_and_swap(p['block2']['JMIN'], p['block2']['JMAX']) + p['block2']['KMIN'], p['block2']['KMAX'] = check_and_swap(p['block2']['KMIN'], p['block2']['KMAX']) + voi = [p['block2']['IMIN'], p['block2']['IMAX'], p['block2']['JMIN'], p['block2']['JMAX'],p['block2']['KMIN'], p['block2']['KMAX']] + CreateSubset(block_source, voi, name='periodic '+str(periodic_indx)) \ No newline at end of file diff --git a/python/examples/translated_periodicity/vspt_tutorial.ipynb b/python/examples/translated_periodicity/vspt_tutorial.ipynb index 4615ac1..de776be 100644 --- a/python/examples/translated_periodicity/vspt_tutorial.ipynb +++ b/python/examples/translated_periodicity/vspt_tutorial.ipynb @@ -89,6 +89,21 @@ "Finding connectivity\n", "gcd to use 16\n" ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Finding nearest blocks comparing 1436 with 1437: 100%|██████████| 1033203/1033203 [51:52<00:00, 331.98it/s]\n", + "Checking connections block 1436 with 1437: 100%|██████████| 8607/8607 [06:04<00:00, 23.63it/s]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Organizing split and outerfaces\n" + ] } ], "source": [ diff --git a/python/plot3d/connectivity.py b/python/plot3d/connectivity.py index b952b37..0532430 100644 --- a/python/plot3d/connectivity.py +++ b/python/plot3d/connectivity.py @@ -193,7 +193,7 @@ def get_face_intersection(face1:Face,face2:Face,block1:Block,block2:Block,tol:fl df = pd.concat([df, pd.DataFrame(match_location)], ignore_index=True) # Checking for split faces - if len(df)>4: + if len(df)>=4: if (__check_edge(df)): df = pd.DataFrame() # If it's an edge else: # not edge @@ -219,7 +219,7 @@ def get_face_intersection(face1:Face,face2:Face,block1:Block,block2:Block,tol:fl df = __filter_block_increasing(df,'j2') # Do a final check after doing all these checks - if len(df)>4: # Greater than 4 because match can occur with simply 4 corners but the interior doesn't match. + if len(df)>=4: # Greater than 4 because match can occur with simply 4 corners but the interior doesn't match. # Check for Split faces ## Block 1 main_face = create_face_from_diagonals(block1,imin=I1[0],imax=I1[1], jmin=J1[0],jmax=J1[1],kmin=K1[0],kmax=K1[1]) diff --git a/python/plot3d/periodicity.py b/python/plot3d/periodicity.py index 04cdfb8..27ad390 100644 --- a/python/plot3d/periodicity.py +++ b/python/plot3d/periodicity.py @@ -570,10 +570,13 @@ def shift_blocks(sign:int=1): pbar = tqdm(total = len(lower_connected_faces)) periodicity_tol = 1E-6 - while len(lower_connected_faces)>0 and times_failed<3: + while len(lower_connected_faces)>0: periodic_found = False face1 = lower_connected_faces[0] - for face2 in upper_connected_faces: + for indx in range(len(upper_connected_faces)): + face2 = upper_connected_faces[indx] + if face1.BlockIndex == 26 and face2.BlockIndex == 62: + print('check') # Check if surfaces are periodic with each other pbar.set_description(f"Checking connections block {face1.blockIndex} with {face2.blockIndex}") # Shift block 1 -> Check periodicity -> if not periodic -> shift Block 1 opposite direction -> Check periodicity @@ -584,8 +587,8 @@ def shift_blocks(sign:int=1): _, periodic_faces_temp, split_faces_temp = __periodicity_check__(face1,face2,block1_shifted, block2,periodicity_tol) if len(periodic_faces_temp) > 0: - lower_connected_faces.remove(face1) - upper_connected_faces.remove(face2) + lower_connected_faces.pop(0) + upper_connected_faces.pop(indx) periodic_faces.append(periodic_faces_temp) periodic_faces_export.append(face_matches_to_dict(periodic_faces_temp[0],periodic_faces_temp[1],block1_shifted,block2)) lower_split_faces = [s for s in split_faces_temp if s.BlockIndex in lower_blocks] @@ -598,9 +601,9 @@ def shift_blocks(sign:int=1): if periodic_found == False: # Lets switch the order - non_matching.append(face1) - lower_connected_faces.remove(face1) - periodicity_tol *=10 + non_matching.append(deepcopy(face1)) + lower_connected_faces.pop(0) + # periodicity_tol *=10 times_failed+=1 if len(non_matching)>0: @@ -743,7 +746,7 @@ def __periodicity_check__(face1:Face, face2:Face,block1:Block,block2:Block,tol:f df,split_face1,split_face2 = get_face_intersection(face1,face2,block1,block2,tol) - if len(df)>4: + if len(df)>=4: f1 = create_face_from_diagonals(block1,imin=df['i1'].min(),jmin=df['j1'].min(),kmin=df['k1'].min(), imax=df['i1'].max(),jmax=df['j1'].max(),kmax=df['k1'].max()) f1.set_block_index(face1.blockIndex) diff --git a/python/pyproject.toml b/python/pyproject.toml index 96c1fac..aba6a60 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "plot3d" -version = "1.5.7" +version = "1.5.8" description = "Plot3D python utilities for reading and writing and also finding connectivity between blocks" authors = ["Paht Juangphanich "]