diff --git a/.gitignore b/.gitignore index 026c945..0c2a4c2 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,5 @@ # Created by venv; see https://docs.python.org/3/library/venv.html ShipD/ -ShipGen/ Scripts/ Lib/ Include/ diff --git a/ShipGen/.gitignore b/ShipGen/.gitignore new file mode 100644 index 0000000..916c6b0 --- /dev/null +++ b/ShipGen/.gitignore @@ -0,0 +1,3 @@ +DEMO_Hulls/ +__pycache__ +data/*.npy \ No newline at end of file diff --git a/ShipGen/C_ShipGen_DEMO.ipynb b/ShipGen/C_ShipGen_DEMO.ipynb new file mode 100644 index 0000000..a94e169 --- /dev/null +++ b/ShipGen/C_ShipGen_DEMO.ipynb @@ -0,0 +1,678 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## C-ShipGen: Sample Tailored Ship Hulls from a Tabular DDPM" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### Set Up Tasks: Don't alter Please #####" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# import the fun\n", + "import sys\n", + "import os \n", + "# sys.path.append('./tools')\n", + "# sys.path.append('./data')\n", + "\n", + "import numpy as np\n", + "from tqdm import tqdm\n", + "import math\n", + "import matplotlib.pyplot as plt\n", + "import torch\n", + "import torch.nn as nn\n", + "import torch.nn.functional as F\n", + "from tools import Guided_Cond_DDPM_Tools as GC_DDPM\n", + "\n", + "from sklearn.decomposition import PCA\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from tools.HullParameterization import Hull_Parameterization as HP\n", + "\n", + "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n", + "\n", + "\n", + "np.set_printoptions(suppress=True) # don't use scientific notation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(82168, 45)\n", + "(82793, 44)\n", + "(44, 2)\n" + ] + } + ], + "source": [ + "# Load in the Data:\n", + "\n", + "#Step 1: Load in the data\n", + "DesVec = np.load('./data/DesVec_82k.npy', allow_pickle=True)\n", + "print(DesVec.shape)\n", + "\n", + "DesVec_neg = np.load('./data/Negative_DesVec_82k.npy', allow_pickle=True)\n", + "print(DesVec_neg.shape)\n", + "\n", + "\n", + "# Now lets clean up X\n", + "\n", + "idx_BBFactors = [33,34,35,36,37]\n", + "idx_BB = 31\n", + "\n", + "idx_SBFactors = [38,39,40,41,42,43,44]\n", + "idx_SB = 32\n", + "\n", + "for i in range(0,len(DesVec)):\n", + " \n", + " DesVec[i,idx_BBFactors] = DesVec[i,idx_BB] * DesVec[i,idx_BBFactors] \n", + " DesVec[i,idx_SBFactors] = DesVec[i,idx_SB] * DesVec[i,idx_SBFactors]\n", + "\n", + "\n", + "\n", + "Y = np.load('./data/GeometricMeasures.npy', allow_pickle=True)\n", + "\n", + "LenRatios = np.load('./data/Length_Ratios.npy', allow_pickle=True)\n", + "\n", + "\n", + "X_LIMITS = np.load('./data/X_LIMITS.npy')\n", + "\n", + "print(X_LIMITS.shape)\n", + "\n", + "X_lower_lim = [X_LIMITS[:,0].tolist()] \n", + "X_upper_lim = [X_LIMITS[:,1].tolist()]\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(array([77257, 77257, 77257, 77257], dtype=int64), array([1, 2, 3, 4], dtype=int64))\n", + "(82168, 101)\n", + "(82168,)\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\joe\\AppData\\Local\\Temp\\ipykernel_4716\\214063027.py:4: RuntimeWarning: invalid value encountered in log10\n", + " VolVec = np.log10(Y[:,1*num_WL_Steps:2*num_WL_Steps])\n" + ] + } + ], + "source": [ + "#Set up Conditioning Vectors:\n", + "num_WL_Steps = 101\n", + "\n", + "VolVec = np.log10(Y[:,1*num_WL_Steps:2*num_WL_Steps])\n", + "idx = np.where(np.isnan(VolVec))\n", + "print(idx)\n", + "\n", + "VolVec[idx] = -6.0 #fix nan to dummy value\n", + "\n", + "print(VolVec.shape)\n", + "\n", + "DdVec = DesVec[:,4]\n", + "BOAVec = np.amax(LenRatios[:,1:3], axis=1)\n", + "print(BOAVec.shape) \n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the file for architecting the network, diffusion parameters, and training\n", + "\n", + "DDPM_Dict = {\n", + " 'xdim' : len(DesVec[0])-1, # Dimension of parametric design vector\n", + " 'datalength': len(DesVec), # number of samples\n", + " 'X_LL' : X_lower_lim, # lower limits of parametric design vector variables\n", + " 'X_UL' : X_upper_lim,\n", + " 'ydim': 0, # Number of objectives\n", + " 'cdim': 4, # number of conditioning inputs\n", + " 'gamma' : 0.2, # weight of feasibility guidance for guided sampling\n", + " 'lambda': [0.3,0.3], # weight of drag guidance for guided sampling\n", + " #'lambdas': [1,1,1,1,1,1,1], # dummy variable for performance guided sampling\n", + " 'tdim': 128, # dimension of latent variable\n", + " 'net': [1024,1024,1024,1024], # network architecture\n", + " 'batch_size': 1024, # batch size\n", + " 'Training_Epochs': 10000, # number of training epochs\n", + " 'Diffusion_Timesteps': 1000, # number of diffusion timesteps\n", + " 'lr' : 0.00025, # learning rate\n", + " 'weight_decay': 0.0, # weight decay\n", + " 'device_name': device} # gpu device name\n", + "\n", + "\n", + "Classify_Dict = {\n", + " 'xdim' : len(DesVec[0])-1,\n", + " 'cdim': 1,\n", + " 'tdim': 128,\n", + " 'net': [64,64,64],\n", + " 'Training_Epochs': 150000,\n", + " 'device_name': device}\n", + "\n", + "nodes = 512\n", + "Drag_Reg_Dict = {\n", + " 'xdim' : len(DesVec[0])-1, # Dimension of parametric design vector\n", + " 'ydim': 1, # trains regression model for each objective\n", + " 'tdim': nodes, # dimension of latent variable\n", + " 'net': [nodes,nodes,nodes,nodes], # network architecture \n", + " 'Training_Epochs': 50000, #30000 # number of training epochs\n", + " 'batch_size': 1024, # batch size\n", + " 'Model_Label': 'Regressor_CT', # labels for regressors \n", + " 'lr' : 0.001, # learning rate\n", + " 'weight_decay': 0.0, # weight decay\n", + " 'device_name': device} \n", + "\n", + "nodes = 256\n", + "LOA_wBulb_Reg_Dict = {\n", + " 'xdim' : len(DesVec[0])-1, # Dimension of parametric design vector\n", + " 'ydim': 1, # trains regression model for each objective\n", + " 'tdim': nodes, # dimension of latent variable\n", + " 'net': [nodes,nodes,nodes], # network architecture \n", + " 'Training_Epochs': 150000, # number of training epochs\n", + " 'batch_size': 1024, # batch size\n", + " 'Model_Label': 'Regressor_LOA_wBulb', # labels for regressors\n", + " \n", + " 'lr' : 0.001, # learning rate\n", + " 'weight_decay': 0.0, # weight decay\n", + " 'device_name': device} \n", + "\n", + "WL_Reg_Dict = {\n", + " \"xdim\": len(DesVec[0]),\n", + " \"ydim\": 1, \n", + " \"tdim\": 512, \n", + " \"net\": [512, 512, 512], \n", + " \"Training_Epochs\": 30000, \n", + " \"batch_size\": 1024, \n", + " \"Model_Label\": \n", + " \"Regressor_WL\", \n", + " \"lr\": 0.001, \n", + " \"weight_decay\": 0.0, \n", + " \"device_name\": device}\n", + "\n", + "Vol_Reg_Dict = {\n", + " \"xdim\": len(DesVec[0]), \n", + " \"ydim\": 1, \n", + " \"tdim\": 512, \n", + " \"net\": [512, 512, 512], \n", + " \"Training_Epochs\": 30000, \n", + " \"batch_size\": 1024, \n", + " \"Model_Label\": \"Regressor_WL\", \n", + " \"lr\": 0.001, \n", + " \"weight_decay\": 0.0, \n", + " \"device_name\": device}\n", + "\n", + "\n", + "\n", + "\n", + "T = GC_DDPM.GuidedDiffusionEnv(DDPM_Dict,\n", + " Classify_Dict,\n", + " Drag_Reg_Dict,\n", + " LOA_wBulb_Reg_Dict,\n", + " WL_Reg_Dict,\n", + " Vol_Reg_Dict,\n", + " X= DesVec[:,1:],\n", + " X_neg= DesVec_neg,\n", + " VolVec = VolVec, \n", + " BOAVec = BOAVec, \n", + " DdVec = DdVec)\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "diffusion_path = './TrainedModels/CShipGen_diffusion.pth'\n", + "T.load_trained_diffusion_model(diffusion_path)\n", + "\n", + "classifier_path = './TrainedModels/Constraint_Classifier_150000Epochs.pth' \n", + "\n", + "T.load_trained_classifier_model(classifier_path)\n", + "\n", + "\n", + "PATHS = ['./TrainedModels/Regressor_CT.pth',\n", + " './TrainedModels/Regressor_LOA_wBulb.pth',\n", + " './TrainedModels/Regressor_WL.pth',\n", + " './TrainedModels/Regressor_Vol.pth']\n", + "\n", + "\n", + "T.load_trained_Drag_regression_models(PATHS)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define the Ship's Principal Characteristics ##" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "#Sample from the Model:\n", + "num_samples = 100\n", + "\n", + "Ship = np.array([333, 42.624, 11.28, 29.064, 97561,16]) #[LOA(m), Beam(m), Draft(m), Depth(m), Volume(m^3), U(m/s)] # This is for an aircraft carrier dimensioned ship\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Generate the Hulls using C-ShipGen ##" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Generating Hulls\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 967/967 [00:41<00:00, 23.20it/s]\n", + "100%|██████████| 32/32 [00:00<00:00, 91.00it/s]" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(100, 45)\n", + "Predicted Mean Drag of Guidance samples: 16247670.0 N\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "\n", + "# Run the Loop on the other samples:\n", + "\n", + "print('Generating Hulls')\n", + "\n", + "LOA = Ship[0] #in meters\n", + "BoL = Ship[1]/LOA #beam to length ratio\n", + "ToD = Ship[2]/Ship[3] #Draft to depth ratio\n", + "DoL = Ship[3]/LOA #Depth to length ratio\n", + "Vol = np.log10(Ship[4]/LOA**3) # to normalize Volume by LOA**3\n", + "\n", + "U = Ship[5] # 12.86 #m/s = 25 knots\n", + "\n", + "dim_d = np.array([[ToD, U, LOA]]) #Drag_conditioning is [ToD, U(m/s), LOA (m)]\n", + "\n", + "drag_cond = np.repeat(dim_d, num_samples, axis=0) #reapeat \n", + "#print(cond.shape)\n", + "dim_g = np.array([[ToD, BoL, DoL, Vol]])\n", + "\n", + "geom_cond = np.repeat(dim_g, num_samples, axis=0) #reapeat \n", + "#print(cond.shape)\n", + "\n", + "X_gen, unnorm = T.gen_vol_drag_guided_samples(geom_cond, drag_cond)\n", + "\n", + "print(X_gen.shape)\n", + "\n", + "\n", + "Rt_guidance = T.Predict_Drag(unnorm, drag_cond)\n", + "Drag_Guidance = np.mean(Rt_guidance)\n", + "\n", + "\n", + "print('Predicted Mean Drag of Guidance samples: ' + str(Drag_Guidance) + ' N')\n", + "#print('Minimum Drag of Guidance samples: ' + str(np.amin(Rt_guidance)) + ' N')\n", + "\n", + "\n", + " \n", + "\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Clean up the Vectors and Check Feasibility ###" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Checking Feasibility of Samples\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 100/100 [00:00<00:00, 1171.57it/s]\n" + ] + }, + { + "data": { + "text/plain": [ + "\"\\nprint(len(valid_idx))\\nsample_vol = []\\nsample_BOA = []\\nsample_Dd = []\\nsample_LOA = []\\nsample_LOA_wBulb = []\\nidx_to_remove = []\\n\\nfor i in tqdm(range(0,len(valid_idx))):\\n hull = HP(x_samples[valid_idx[i]]) \\n #print(i)\\n try:\\n Z = hull.Calc_VolumeProperties(NUM_WL = 101, PointsPerWL = 1000) \\n sample_vol.append(HP.interp(hull.Volumes, Z, Ship[2])) #interpolate to the draft of the ship\\n BOA = max(hull.Calc_Max_Beam_midship(), hull.Calc_Max_Beam_PC())\\n sample_BOA.append(BOA)\\n sample_Dd.append(hull.Dd)\\n sample_LOA.append(hull.LOA)\\n sample_LOA_wBulb.append(hull.Calc_LOA_wBulb())\\n except:\\n print('error at hull {}'.format(i))\\n idx_to_remove.append(i)\\n\\n continue\\n\\n#Remove the samples that failed to calculate volume properties\\nvalid_idx = np.delete(valid_idx, idx_to_remove)\\nprint(len(valid_idx))\\n\"" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\n", + "x_samples = X_gen\n", + "\n", + "#print(x_samples[0:3])\n", + " \n", + "print('Checking Feasibility of Samples')\n", + "\n", + "for i in range(0,len(x_samples)):\n", + " \n", + " x_samples[i,idx_BB] = (x_samples[i,idx_BB] + 0.5) // 1 #int rounds to 1 or 0\n", + " x_samples[i,idx_SB] = (x_samples[i,idx_SB] + 0.5) // 1 #int rounds to 1 or 0\n", + " \n", + " \n", + " x_samples[i,idx_BBFactors] = x_samples[i,idx_BB] * x_samples[i,idx_BBFactors] \n", + " x_samples[i,idx_SBFactors] = x_samples[i,idx_SB] * x_samples[i,idx_SBFactors]\n", + "\n", + "\n", + "\n", + "#Check the constraint violations for the sampled designs\n", + "constraints = []\n", + "sum_violation = []\n", + "cons = []\n", + "valid_idx = []\n", + "\n", + "for i in tqdm(range(0,len(x_samples))):\n", + " hull = HP(x_samples[i])\n", + " constraints.append(hull.input_Constraints())\n", + " cons.append(constraints[i] > 0)\n", + " if sum(cons[i]) == 0:\n", + " valid_idx.append(i)\n", + " #hull.Calc_VolumeProperties(NUM_WL = 101, PointsPerWL = 1000)\n", + " sum_violation.append(sum(cons[i]))\n", + "\n", + "'''\n", + "print(len(valid_idx))\n", + "sample_vol = []\n", + "sample_BOA = []\n", + "sample_Dd = []\n", + "sample_LOA = []\n", + "sample_LOA_wBulb = []\n", + "idx_to_remove = []\n", + "\n", + "for i in tqdm(range(0,len(valid_idx))):\n", + " hull = HP(x_samples[valid_idx[i]]) \n", + " #print(i)\n", + " try:\n", + " Z = hull.Calc_VolumeProperties(NUM_WL = 101, PointsPerWL = 1000) \n", + " sample_vol.append(HP.interp(hull.Volumes, Z, Ship[2])) #interpolate to the draft of the ship\n", + " BOA = max(hull.Calc_Max_Beam_midship(), hull.Calc_Max_Beam_PC())\n", + " sample_BOA.append(BOA)\n", + " sample_Dd.append(hull.Dd)\n", + " sample_LOA.append(hull.LOA)\n", + " sample_LOA_wBulb.append(hull.Calc_LOA_wBulb())\n", + " except:\n", + " print('error at hull {}'.format(i))\n", + " idx_to_remove.append(i)\n", + "\n", + " continue\n", + "\n", + "#Remove the samples that failed to calculate volume properties\n", + "valid_idx = np.delete(valid_idx, idx_to_remove)\n", + "print(len(valid_idx))\n", + "'''\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Generate Stl of valid Hull Designs ###\n", + "\n", + "Note: \n", + "\n", + "Not all generated samples are feasible since C-ShipGen is a statistical Model.\n", + "\n", + "Similarly, C-ShipGen does not generate hull designs exactly to the dimensions specified by the user; however, these designs are close to the intended dimensions. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Designs with Minimum Drag: \n", + "Demo_Hull_29\n", + "Demo_Hull_73\n", + "Demo_Hull_15\n", + "Demo_Hull_41\n", + "Demo_Hull_62\n", + "Demo_Hull_26\n", + "Demo_Hull_57\n", + "Demo_Hull_83\n", + "Demo_Hull_1\n", + "Demo_Hull_47\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 10/10 [00:02<00:00, 4.50it/s]" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Minimum Predicted Drag: 6978521.5 N\n", + "Design with Minimum Predicted Drag: \n", + "Demo_Hull_29\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "path = \"./DEMO_Hulls/\"\n", + "\n", + "label = 'Demo_Hull_'\n", + "\n", + "isExist = os.path.exists(path)\n", + "if not isExist:\n", + " os.makedirs(path)\n", + "\n", + "# print the indices of the 10 hulls with the lowest drag\n", + "print('Designs with Minimum Drag: ')\n", + "idxs = np.argsort(Rt_guidance[valid_idx].flatten())\n", + "for i in range(10):\n", + " print(label + str(valid_idx[idxs[i]]))\n", + "\n", + "for i in tqdm(range(0,10)):\n", + " Hull = HP(x_samples[valid_idx[idxs[i]]])\n", + " #make the .stl file of the hull:\n", + " strpath = path+label + '_' + str(valid_idx[idxs[i]])\n", + " try:\n", + " mesh = Hull.gen_stl(NUM_WL=47, PointsPerWL=151, bit_AddTransom = 1, bit_AddDeckLid = 1, bit_RefineBowAndStern = 1,namepath = strpath)\n", + " except:\n", + " print('Error at hull {}'.format(valid_idx[idxs[i]]))\n", + "\n", + "idx_min_pred_drag = np.argmin(Rt_guidance[valid_idx])\n", + "print('Minimum Predicted Drag: ' + str(Rt_guidance[valid_idx][idx_min_pred_drag][0]) + ' N')\n", + "print('Design with Minimum Predicted Drag: ')\n", + "print(label + str(valid_idx[idx_min_pred_drag])) #Highlight design with minimum predicted drag\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\"\\nLabels = ['Low', 'Medium', 'High']\\nWL = [47, 111, 203]\\nPPW = [151, 301, 601]\\n\\nidx = 87\\n\\nfor i in tqdm(range(len(Labels))):\\n Hull = HP(x_samples[idx])\\n #make the .stl file of the hull:\\n strpath = path+label + '_' + str(idx) + '_' + Labels[i]\\n try:\\n mesh = Hull.gen_stl(NUM_WL=WL[i], PointsPerWL=PPW[i], bit_AddTransom = 1, bit_AddDeckLid = 1, bit_RefineBowAndStern = 1,namepath = strpath)\\n except:\\n print('Error at hull {}'.format(valid_idx[idxs[i]]))\\n\"" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Design a hull with multiple fidelity levels\n", + "'''\n", + "Labels = ['Low', 'Medium', 'High']\n", + "WL = [47, 111, 203]\n", + "PPW = [151, 301, 601]\n", + "\n", + "idx = 87\n", + "\n", + "for i in tqdm(range(len(Labels))):\n", + " Hull = HP(x_samples[idx])\n", + " #make the .stl file of the hull:\n", + " strpath = path+label + '_' + str(idx) + '_' + Labels[i]\n", + " try:\n", + " mesh = Hull.gen_stl(NUM_WL=WL[i], PointsPerWL=PPW[i], bit_AddTransom = 1, bit_AddDeckLid = 1, bit_RefineBowAndStern = 1,namepath = strpath)\n", + " except:\n", + " print('Error at hull {}'.format(valid_idx[idxs[i]]))\n", + "'''" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "45\n", + "Predicted Drag of Test Hull: 12658716.0 N\n" + ] + } + ], + "source": [ + "# Compare to Test Hull:\n", + "\n", + "design = [[333,\t0.587059958,\t0.409996119,\t0.122972918,\t0.085,\t0.71417088,\t0.211,\t0.120262544,\t10.70168725,\t0.264897988,\t0.961214759,\t-0.27,\t0.15,\t0.01,\t0.3,\t2.172250534,\t-2.220618947,\t0,\t0.1,\t0.05,\t0,\t0,\t0.104249968,\t0.37043854,\t0.006300723,\t-2.478684046,\t2.882864263,\t3.320624286,\t0.076330577,\t0.48538725,\t0.455186006,\t1,\t1,\t0.01,\t0.4,\t0.17155627,\t0.38035235,\t0.269331342,\t0.7,\t0.01,\t1,\t0.7,\t0.025,\t0.99,\t0.147764723]]\n", + "\n", + "print(len(design[0]))\n", + "\n", + "hull = HP(design[0])\n", + "\n", + "strpath = path+'Test_Hull_Nimitz'\n", + "mesh = hull.gen_stl(NUM_WL=47, PointsPerWL=151, bit_AddTransom = 1, bit_AddDeckLid = 1, bit_RefineBowAndStern = 1,namepath = strpath)\n", + "\n", + "unnormed_des = T.data_norm.transform_Data(design[0][1:])\n", + "#unnormed_des = torch.from_numpy(unnormed_des.astype('float32'))\n", + "\n", + "Rt_Test = T.Predict_Drag(unnormed_des, drag_cond[0:1])\n", + "\n", + "print('Predicted Drag of Test Hull: ' + str(Rt_Test[0,0]) + ' N')\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/ShipGen/Generate.py b/ShipGen/Generate.py new file mode 100644 index 0000000..321ce90 --- /dev/null +++ b/ShipGen/Generate.py @@ -0,0 +1,260 @@ +import os + +import numpy as np +from tqdm import tqdm +import torch +from .tools import Guided_Cond_DDPM_Tools as GC_DDPM +from .tools.HullParameterization import Hull_Parameterization as HP +from stl import Mesh + +data_path = os.path.normpath(os.path.join(__file__,'../data')) +trained_models_path = os.path.normpath(os.path.join(__file__,'../TrainedModels')) + +device = torch.device("cuda" if torch.cuda.is_available() else "cpu") + +class ShipCharacteristics: + LOA:float + Beam:float + Draft:float + Depth:float + Volume:float + U:float + def __init__(self, + LOA:float, + Beam:float, + Draft:float, + Depth:float, + Volume:float, + U:float): + self.LOA = LOA + self.Beam = Beam + self.Draft = Draft + self.Depth = Depth + self.Volume = Volume + self.U = U + def as_array(self) -> np.ndarray: + return np.array([self.LOA, self.Beam, self.Draft, self.Depth, self.Volume, self.U]) + +def generate_hulls(find_best:int, + out_of:int|None = None, + characteristics:ShipCharacteristics|None = None, + out_of_ratio:int = 10) -> list[Mesh]: + if characteristics is None: + Ship = np.array([333, 42.624, 11.28, 29.064, 97561,16]) + else: Ship = characteristics.as_array() + if out_of is None: num_samples = find_best * out_of_ratio + else: num_samples = out_of + + LOA:float = Ship[0] #in meters + BoL:float = Ship[1]/LOA #beam to length ratio + ToD:float = Ship[2]/Ship[3] #Draft to depth ratio + DoL:float = Ship[3]/LOA #Depth to length ratio + Vol:float = np.log10(Ship[4]/LOA**3) # to normalize Volume by LOA**3 + + U:float = Ship[5] # 12.86 #m/s = 25 knots + + dim_d = np.array([[ToD, U, LOA]]) #Drag_conditioning is [ToD, U(m/s), LOA (m)] + + drag_cond = np.repeat(dim_d, num_samples, axis=0) #reapeat + dim_g = np.array([[ToD, BoL, DoL, Vol]]) + + geom_cond = np.repeat(dim_g, num_samples, axis=0) #reapeat + + x_samples, unnorm = T.gen_vol_drag_guided_samples(geom_cond, drag_cond) + + Rt_guidance = T.Predict_Drag(unnorm, drag_cond) + + valid_idx = check_feasibility(x_samples) + + idxs = np.argsort(Rt_guidance[valid_idx].flatten()) + meshs:list[Mesh] = [] + for i in tqdm(range(0,find_best)): + Hull = HP(x_samples[valid_idx[idxs[i]]]) + #make the .stl file of the hull: + try: + meshs.append(Hull.gen_stl(NUM_WL=47, PointsPerWL=151, bit_AddTransom = 1, bit_AddDeckLid = 1, bit_RefineBowAndStern = 1,namepath = None)) + except: + print('Error at hull {}'.format(valid_idx[idxs[i]])) + return meshs + +def mesh_to_triangles(mesh:Mesh) -> np.ndarray: + return np.array([mesh.v0, mesh.v1, mesh.v2]).swapaxes(0, 1) +def mesh_list_to_triangles(mesh_list:list[Mesh]) -> list[np.ndarray]: + return list(map(mesh_to_triangles, mesh_list)) + +def normalize_x(mesh:Mesh, scale:float = 10): + mesh.vectors -= np.array([mesh.min_[0], 0, 0]) + mesh.vectors *= scale / (mesh.max_[0] - mesh.min_[0]) + +def load_data(): + # Load in the Data: + + #Step 1: Load in the data + DesVec = np.load(f'{data_path}/DesVec_82k.npy', allow_pickle=True) + + DesVec_neg = np.load(f'{data_path}/Negative_DesVec_82k.npy', allow_pickle=True) + + # Now lets clean up X + + idx_BBFactors = [33,34,35,36,37] + idx_BB = 31 + + idx_SBFactors = [38,39,40,41,42,43,44] + idx_SB = 32 + + for i in range(0,len(DesVec)): + DesVec[i,idx_BBFactors] = DesVec[i,idx_BB] * DesVec[i,idx_BBFactors] + DesVec[i,idx_SBFactors] = DesVec[i,idx_SB] * DesVec[i,idx_SBFactors] + + + Y = np.load(f'{data_path}/GeometricMeasures.npy', allow_pickle=True) + + LenRatios = np.load(f'{data_path}/Length_Ratios.npy', allow_pickle=True) + + X_LIMITS = np.load(f'{data_path}/X_LIMITS.npy') + + #Set up Conditioning Vectors: + num_WL_Steps = 101 + + old_err = np.seterr(all='ignore') + VolVec = np.log10(Y[:,1*num_WL_Steps:2*num_WL_Steps]) + idx = np.where(np.isnan(VolVec)) + VolVec[idx] = -6.0 #fix nan to dummy value + np.seterr(**old_err) + + + BOAVec = np.amax(LenRatios[:,1:3], axis=1) + return ( + (idx_BB, idx_SB, idx_BBFactors, idx_SBFactors), + (DesVec, X_LIMITS, DesVec_neg, VolVec, BOAVec) + ) + +def create_network(DesVec, X_LIMITS, DesVec_neg, VolVec, BOAVec) -> GC_DDPM.GuidedDiffusionEnv: + # Set up the file for architecting the network, diffusion parameters, and training + X_lower_lim = [X_LIMITS[:,0].tolist()] + X_upper_lim = [X_LIMITS[:,1].tolist()] + DDPM_Dict = { + 'xdim' : len(DesVec[0])-1, # Dimension of parametric design vector + 'datalength': len(DesVec), # number of samples + 'X_LL' : X_lower_lim, # lower limits of parametric design vector variables + 'X_UL' : X_upper_lim, + 'ydim': 0, # Number of objectives + 'cdim': 4, # number of conditioning inputs + 'gamma' : 0.2, # weight of feasibility guidance for guided sampling + 'lambda': [0.3,0.3], # weight of drag guidance for guided sampling + #'lambdas': [1,1,1,1,1,1,1], # dummy variable for performance guided sampling + 'tdim': 128, # dimension of latent variable + 'net': [1024,1024,1024,1024], # network architecture + 'batch_size': 1024, # batch size + 'Training_Epochs': 10000, # number of training epochs + 'Diffusion_Timesteps': 1000, # number of diffusion timesteps + 'lr' : 0.00025, # learning rate + 'weight_decay': 0.0, # weight decay + 'device_name': device} # gpu device name + + Classify_Dict = { + 'xdim' : len(DesVec[0])-1, + 'cdim': 1, + 'tdim': 128, + 'net': [64,64,64], + 'Training_Epochs': 150000, + 'device_name': device} + + nodes = 512 + Drag_Reg_Dict = { + 'xdim' : len(DesVec[0])-1, # Dimension of parametric design vector + 'ydim': 1, # trains regression model for each objective + 'tdim': nodes, # dimension of latent variable + 'net': [nodes,nodes,nodes,nodes], # network architecture + 'Training_Epochs': 50000, #30000 # number of training epochs + 'batch_size': 1024, # batch size + 'Model_Label': 'Regressor_CT', # labels for regressors + 'lr' : 0.001, # learning rate + 'weight_decay': 0.0, # weight decay + 'device_name': device} + + nodes = 256 + LOA_wBulb_Reg_Dict = { + 'xdim' : len(DesVec[0])-1, # Dimension of parametric design vector + 'ydim': 1, # trains regression model for each objective + 'tdim': nodes, # dimension of latent variable + 'net': [nodes,nodes,nodes], # network architecture + 'Training_Epochs': 150000, # number of training epochs + 'batch_size': 1024, # batch size + 'Model_Label': 'Regressor_LOA_wBulb', # labels for regressors + + 'lr' : 0.001, # learning rate + 'weight_decay': 0.0, # weight decay + 'device_name': device} + + WL_Reg_Dict = { + "xdim": len(DesVec[0]), + "ydim": 1, + "tdim": 512, + "net": [512, 512, 512], + "Training_Epochs": 30000, + "batch_size": 1024, + "Model_Label": + "Regressor_WL", + "lr": 0.001, + "weight_decay": 0.0, + "device_name": device} + + Vol_Reg_Dict = { + "xdim": len(DesVec[0]), + "ydim": 1, + "tdim": 512, + "net": [512, 512, 512], + "Training_Epochs": 30000, + "batch_size": 1024, + "Model_Label": "Regressor_WL", + "lr": 0.001, + "weight_decay": 0.0, + "device_name": device} + + T = GC_DDPM.GuidedDiffusionEnv(DDPM_Dict, + Classify_Dict, + Drag_Reg_Dict, + LOA_wBulb_Reg_Dict, + WL_Reg_Dict, + Vol_Reg_Dict, + X= DesVec[:,1:], + X_neg= DesVec_neg, + VolVec = VolVec, + BOAVec = BOAVec, + DdVec = DesVec[:,4]) + + diffusion_path = f'{trained_models_path}/CShipGen_diffusion.pth' + T.load_trained_diffusion_model(diffusion_path) + classifier_path = f'{trained_models_path}/Constraint_Classifier_150000Epochs.pth' + T.load_trained_classifier_model(classifier_path) + + PATHS = [f'{trained_models_path}/Regressor_CT.pth', + f'{trained_models_path}/Regressor_LOA_wBulb.pth', + f'{trained_models_path}/Regressor_WL.pth', + f'{trained_models_path}/Regressor_Vol.pth'] + + T.load_trained_Drag_regression_models(PATHS) + return T +feasibility_data, network_data = load_data() +T = create_network(*network_data) + +def check_feasibility(x_samples:np.ndarray): + idx_BB, idx_SB, idx_BBFactors, idx_SBFactors = feasibility_data + for i in range(0,len(x_samples)): + x_samples[i,idx_BB] = (x_samples[i,idx_BB] + 0.5) // 1 #int rounds to 1 or 0 + x_samples[i,idx_SB] = (x_samples[i,idx_SB] + 0.5) // 1 #int rounds to 1 or 0 + + + x_samples[i,idx_BBFactors] = x_samples[i,idx_BB] * x_samples[i,idx_BBFactors] + x_samples[i,idx_SBFactors] = x_samples[i,idx_SB] * x_samples[i,idx_SBFactors] + + #Check the constraint violations for the sampled designs + valid_idx:list[int] = [] + + for i in tqdm(range(0,len(x_samples))): + hull = HP(x_samples[i]) + if sum(hull.input_Constraints() > 0) == 0: + valid_idx.append(i) + return valid_idx + diff --git a/ShipGen/LICENSE b/ShipGen/LICENSE new file mode 100644 index 0000000..2684cf4 --- /dev/null +++ b/ShipGen/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2024 Noah J. Bagazinski + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/ShipGen/TrainedModels/CShipGen_Test.json b/ShipGen/TrainedModels/CShipGen_Test.json new file mode 100644 index 0000000..6fed00f --- /dev/null +++ b/ShipGen/TrainedModels/CShipGen_Test.json @@ -0,0 +1 @@ +{"xdim": 44, "datalength": 82168, "X_LL": [[0.05, 0.0, 0.08333, 0.05, 0.0, 0.05, 0.05, 0.0, 0.0, -1.0, -4.0, -4.0, 0.0, 0.0, -4.0, -4.0, -4.0, -4.0, 0.0, 0.0, 0.0, -3.0, 0.0, 0.0, -4.0, -4.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0]], "X_UL": [[0.9, 0.9, 0.333, 0.25, 1.0, 0.8, 0.5, 45.0, 1.0, 1.0, 4.0, 4.0, 1.0, 1.0, 4.0, 4.0, 4.0, 4.0, 60.0, 1.0, 1.0, 5.0, 1.0, 1.0, 4.0, 4.0, 60.0, 0.5, 0.5, 0.5, 1.0, 1.0, 0.2, 1.0, 1.0, 1.0, 0.33, 1.0, 0.2, 1.0, 1.0, 1.0, 1.0, 0.33]], "ydim": 0, "cdim": 4, "gamma": 0.5, "tdim": 128, "net": [1024, 1024, 1024, 1024], "batch_size": 1024, "Training_Epochs": 10000, "Diffusion_Timesteps": 1000, "lr": 0.00025, "weight_decay": 0.0, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/TrainedModels/CShipGen_Test_diffusion.pth b/ShipGen/TrainedModels/CShipGen_Test_diffusion.pth new file mode 100644 index 0000000..b05cd02 Binary files /dev/null and b/ShipGen/TrainedModels/CShipGen_Test_diffusion.pth differ diff --git a/ShipGen/TrainedModels/CShipGen_diffusion.pth b/ShipGen/TrainedModels/CShipGen_diffusion.pth new file mode 100644 index 0000000..38d48a2 Binary files /dev/null and b/ShipGen/TrainedModels/CShipGen_diffusion.pth differ diff --git a/ShipGen/TrainedModels/CShipGen_diffusion_Dict.json b/ShipGen/TrainedModels/CShipGen_diffusion_Dict.json new file mode 100644 index 0000000..2d67e44 --- /dev/null +++ b/ShipGen/TrainedModels/CShipGen_diffusion_Dict.json @@ -0,0 +1 @@ +{"xdim": 44, "datalength": 82168, "X_LL": [[0.05, 0.0, 0.08333, 0.05, 0.0, 0.05, 0.05, 0.0, 0.0, -1.0, -4.0, -4.0, 0.0, 0.0, -4.0, -4.0, -4.0, -4.0, 0.0, 0.0, 0.0, -3.0, 0.0, 0.0, -4.0, -4.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0]], "X_UL": [[0.9, 0.9, 0.333, 0.25, 1.0, 0.8, 0.5, 45.0, 1.0, 1.0, 4.0, 4.0, 1.0, 1.0, 4.0, 4.0, 4.0, 4.0, 60.0, 1.0, 1.0, 5.0, 1.0, 1.0, 4.0, 4.0, 60.0, 0.5, 0.5, 0.5, 1.0, 1.0, 0.2, 1.0, 1.0, 1.0, 0.33, 1.0, 0.2, 1.0, 1.0, 1.0, 1.0, 0.33]], "ydim": 0, "cdim": 4, "tdim": 128, "net": [1024, 1024, 1024, 1024], "batch_size": 1024, "Training_Epochs": 100000, "Diffusion_Timesteps": 1000, "lr": 0.00025, "weight_decay": 0.0, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/TrainedModels/Constraint_Classifier_150000Epochs.json b/ShipGen/TrainedModels/Constraint_Classifier_150000Epochs.json new file mode 100644 index 0000000..a0015f0 --- /dev/null +++ b/ShipGen/TrainedModels/Constraint_Classifier_150000Epochs.json @@ -0,0 +1 @@ +{"xdim": 44, "cdim": 1, "tdim": 128, "net": [64, 64, 64], "Training_Epochs": 150000, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/TrainedModels/Constraint_Classifier_150000Epochs.pth b/ShipGen/TrainedModels/Constraint_Classifier_150000Epochs.pth new file mode 100644 index 0000000..b8ff6d9 Binary files /dev/null and b/ShipGen/TrainedModels/Constraint_Classifier_150000Epochs.pth differ diff --git a/ShipGen/TrainedModels/Regressor_BOA.pth b/ShipGen/TrainedModels/Regressor_BOA.pth new file mode 100644 index 0000000..6275679 Binary files /dev/null and b/ShipGen/TrainedModels/Regressor_BOA.pth differ diff --git a/ShipGen/TrainedModels/Regressor_BOA_Dict.json b/ShipGen/TrainedModels/Regressor_BOA_Dict.json new file mode 100644 index 0000000..89b2dcd --- /dev/null +++ b/ShipGen/TrainedModels/Regressor_BOA_Dict.json @@ -0,0 +1 @@ +{"xdim": 44, "ydim": 1, "tdim": 256, "net": [256, 256, 256], "Training_Epochs": 1500, "batch_size": 1024, "Model_Label": "Regressor_BOA", "lr": 0.001, "weight_decay": 0.0, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/TrainedModels/Regressor_CT.pth b/ShipGen/TrainedModels/Regressor_CT.pth new file mode 100644 index 0000000..b2c559b Binary files /dev/null and b/ShipGen/TrainedModels/Regressor_CT.pth differ diff --git a/ShipGen/TrainedModels/Regressor_CT_Dict.json b/ShipGen/TrainedModels/Regressor_CT_Dict.json new file mode 100644 index 0000000..82ef6e8 --- /dev/null +++ b/ShipGen/TrainedModels/Regressor_CT_Dict.json @@ -0,0 +1 @@ +{"xdim": 44, "ydim": 1, "tdim": 512, "net": [512, 512, 512, 512], "Training_Epochs": 50000, "batch_size": 1024, "Model_Label": "Regressor_CT", "lr": 0.001, "weight_decay": 0.0, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/TrainedModels/Regressor_LOA_wBulb.pth b/ShipGen/TrainedModels/Regressor_LOA_wBulb.pth new file mode 100644 index 0000000..600fb26 Binary files /dev/null and b/ShipGen/TrainedModels/Regressor_LOA_wBulb.pth differ diff --git a/ShipGen/TrainedModels/Regressor_LOA_wBulb_Dict.json b/ShipGen/TrainedModels/Regressor_LOA_wBulb_Dict.json new file mode 100644 index 0000000..7745a05 --- /dev/null +++ b/ShipGen/TrainedModels/Regressor_LOA_wBulb_Dict.json @@ -0,0 +1 @@ +{"xdim": 44, "ydim": 1, "tdim": 256, "net": [256, 256, 256], "Training_Epochs": 150000, "batch_size": 1024, "Model_Label": "Regressor_LOA_wBulb", "lr": 0.001, "weight_decay": 0.0, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/TrainedModels/Regressor_Vol.pth b/ShipGen/TrainedModels/Regressor_Vol.pth new file mode 100644 index 0000000..1c57b6a Binary files /dev/null and b/ShipGen/TrainedModels/Regressor_Vol.pth differ diff --git a/ShipGen/TrainedModels/Regressor_Vol_Dict.json b/ShipGen/TrainedModels/Regressor_Vol_Dict.json new file mode 100644 index 0000000..91d5509 --- /dev/null +++ b/ShipGen/TrainedModels/Regressor_Vol_Dict.json @@ -0,0 +1 @@ +{"xdim": 44, "ydim": 1, "tdim": 512, "net": [512, 512, 512], "Training_Epochs": 30000, "batch_size": 1024, "Model_Label": "Regressor_Vol", "lr": 0.001, "weight_decay": 0.0, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/TrainedModels/Regressor_WL.pth b/ShipGen/TrainedModels/Regressor_WL.pth new file mode 100644 index 0000000..f69cb4d Binary files /dev/null and b/ShipGen/TrainedModels/Regressor_WL.pth differ diff --git a/ShipGen/TrainedModels/Regressor_WL_Dict.json b/ShipGen/TrainedModels/Regressor_WL_Dict.json new file mode 100644 index 0000000..7c3d8f6 --- /dev/null +++ b/ShipGen/TrainedModels/Regressor_WL_Dict.json @@ -0,0 +1 @@ +{"xdim": 44, "ydim": 1, "tdim": 512, "net": [512, 512, 512], "Training_Epochs": 30000, "batch_size": 1024, "Model_Label": "Regressor_WL", "lr": 0.001, "weight_decay": 0.0, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/TrainedModels/Regressor_WSA.pth b/ShipGen/TrainedModels/Regressor_WSA.pth new file mode 100644 index 0000000..280f6ad Binary files /dev/null and b/ShipGen/TrainedModels/Regressor_WSA.pth differ diff --git a/ShipGen/TrainedModels/Regressor_WSA_Dict.json b/ShipGen/TrainedModels/Regressor_WSA_Dict.json new file mode 100644 index 0000000..c66b79d --- /dev/null +++ b/ShipGen/TrainedModels/Regressor_WSA_Dict.json @@ -0,0 +1 @@ +{"xdim": 44, "ydim": 1, "tdim": 512, "net": [512, 512, 512], "Training_Epochs": 30000, "batch_size": 1024, "Model_Label": "Regressor_WSA", "lr": 0.001, "weight_decay": 0.0, "device_name": "cuda:0"} \ No newline at end of file diff --git a/ShipGen/__init__.py b/ShipGen/__init__.py new file mode 100644 index 0000000..1d07d5d --- /dev/null +++ b/ShipGen/__init__.py @@ -0,0 +1,7 @@ +from .Generate import ( + generate_hulls, + ShipCharacteristics, + mesh_to_triangles, + mesh_list_to_triangles, + normalize_x, +) \ No newline at end of file diff --git a/ShipGen/data/README.txt b/ShipGen/data/README.txt new file mode 100644 index 0000000..8235f6b --- /dev/null +++ b/ShipGen/data/README.txt @@ -0,0 +1 @@ +Please download the relevant files from https://www.dropbox.com/scl/fo/f6bc1kvtlfzped81e2gqu/AJQsD8Ye7kJvkIzi8PTqBdM/C-ShipGen/data?rlkey=akdtdkem92fxqhduzpp85yh95 \ No newline at end of file diff --git a/ShipGen/tools/Guided_Cond_DDPM_Tools.py b/ShipGen/tools/Guided_Cond_DDPM_Tools.py new file mode 100644 index 0000000..771f96b --- /dev/null +++ b/ShipGen/tools/Guided_Cond_DDPM_Tools.py @@ -0,0 +1,1270 @@ +# This script provides a set of tools for creating a guided and/or conditional tabular DDPM Model: + +import numpy as np +import json +import math +import matplotlib.pyplot as plt +import torch +import torch.nn as nn +import torch.nn.functional as F +from tqdm import tqdm +from sklearn.metrics import f1_score +from sklearn.metrics import r2_score + +import sklearn.preprocessing as PP + +''' +========================================== +Set up the data normalizer class +========================================== + +''' + +class Data_Normalizer: + def __init__(self, X_LL_Scaled, X_UL_Scaled,datalength): + + self.normalizer = PP.QuantileTransformer( + output_distribution='normal', + n_quantiles=max(min(datalength // 30, 1000), 10), + subsample=int(1e9) + ) + + self.X_LL_Scaled = X_LL_Scaled + self.X_UL_Scaled = X_UL_Scaled + + self.X_LL_norm = np.zeros((1,len(X_LL_Scaled))) + self.X_UL_norm = np.zeros((1,len(X_LL_Scaled))) + + self.X_mean = np.zeros((1,len(X_LL_Scaled))) + self.X_std = np.zeros((1,len(X_LL_Scaled))) + + def fit_Data(self,X): + + + + x = 2.0*(X-self.X_LL_Scaled)/(self.X_UL_Scaled- self.X_LL_Scaled) - 1.0 + + self.normalizer.fit(x) + x = self.normalizer.transform(x) # Scale Dataset between + #x = (X-self.X_LL_Scaled)/(self.X_UL_Scaled- self.X_LL_Scaled) + + + return x + + def transform_Data(self,X): + x = 2.0*(X-self.X_LL_Scaled)/(self.X_UL_Scaled- self.X_LL_Scaled) - 1.0 + + + x = self.normalizer.transform(x) + return x + + + def scale_X(self,z): + #rescales data + z = self.normalizer.inverse_transform(z) + scaled = (z + 1.0) * 0.5 * (self.X_UL_Scaled - self.X_LL_Scaled) + self.X_LL_Scaled + #scaled = z* (self.X_UL_Scaled - self.X_LL_Scaled) + self.X_LL_Scaled + + ''' + x = self.normalizer.inverse_transform(x) + + #scaled = x* (self.X_UL_norm - self.X_LL_norm) + self.X_LL_norm + ''' + #z = (z + 1.0) * 0.5 * (8.0) + 4.0 + + #scaled = z*self.X_std + self.X_mean + #scaled = self.normalizer.inverse_transform(scaled) + return scaled + +''' +================================================================= +Classifier and Regression Classes +================================================================= +''' +# First Step: make a classifier object: +class Classifier_Model(torch.nn.Module): + def __init__(self, Dict): + nn.Module.__init__(self) + + self.xdim = Dict['xdim'] + self.tdim = Dict['tdim'] + self.cdim = Dict['cdim'] + + self.net = Dict['net'] + self.epochs = Dict['Training_Epochs'] + + self.fc = nn.ModuleList() + + + self.time_embed = nn.Sequential( + nn.Linear(self.tdim, self.tdim), + nn.SiLU(), + nn.Linear(self.tdim, self.tdim)) + + self.X_embed = nn.Linear(self.xdim, self.tdim) + + self.fc.append(self.LinLayer(self.tdim,self.net[0])) + ''' + self.fc.append(self.LinLayer(self.xdim,self.net[0])) + ''' + + for i in range(1, len(self.net)): + self.fc.append(self.LinLayer(self.net[i-1],self.net[i])) + + + self.fc.append(nn.Sequential(nn.Linear(self.net[-1], self.cdim), nn.Sigmoid())) + + def LinLayer(self, dimi, dimo): + + return nn.Sequential(nn.Linear(dimi,dimo), + nn.SiLU(), + #nn.BatchNorm1d(dimo), + nn.Dropout(p=0.1)) + + + def forward(self, x): + + x = self.X_embed(x) + + for i in range(0,len(self.fc)): + x = self.fc[i](x) + + + return x + +class Regression_ResNet(torch.nn.Module): + def __init__(self, Reg_Dict): + nn.Module.__init__(self) + + self.xdim = Reg_Dict['xdim'] + self.ydim = 1 + self.tdim = Reg_Dict['tdim'] + self.net = Reg_Dict['net'] + + self.fc = nn.ModuleList() + + self.fc.append(self.LinLayer(self.tdim,self.net[0])) + + for i in range(1, len(self.net)): + self.fc.append(self.LinLayer(self.net[i-1],self.net[i])) + + self.fc.append(self.LinLayer(self.net[-1], self.tdim)) + ''' + #self.tc = nn.ModuleList() + + #for i in range(0, len(self.net)): + self.tc.append(self.LinLayer(self.tdim,self.net[i])) + self.tc.append(self.LinLayer(self.tdim, self.tdim)) + ''' + self.finalLayer = nn.Sequential(nn.Linear(self.tdim, self.ydim)) + + + self.X_embed = nn.Linear(self.xdim, self.tdim) + #self.T_embed = nn.Linear(self.ydim, self.tdim) + + + def LinLayer(self, dimi, dimo): + + return nn.Sequential(nn.Linear(dimi,dimo), + nn.SiLU(), + nn.LayerNorm(dimo), + nn.Dropout(p=0.1)) + + def forward(self, x): + x = self.X_embed(x) + + res_x = x + + for i in range(0,len(self.fc)): + x = self.fc[i](x) + + x = torch.add(x,res_x) + x = self.finalLayer(x) + + return x + + +class Drag_Regression_ResNet(torch.nn.Module): + def __init__(self, Reg_Dict): + nn.Module.__init__(self) + + self.xdim = Reg_Dict['xdim']+3 # Add 3 Draft, Velocity (Froude Number), and Length scale (LOA) + self.ydim = 1 + self.tdim = Reg_Dict['tdim'] + self.net = Reg_Dict['net'] + + self.fc = nn.ModuleList() + + self.fc.append(self.LinLayer(self.tdim,self.net[0])) + + for i in range(1, len(self.net)): + self.fc.append(self.LinLayer(self.net[i-1],self.net[i])) + + self.fc.append(self.LinLayer(self.net[-1], self.tdim)) + ''' + #self.tc = nn.ModuleList() + + #for i in range(0, len(self.net)): + self.tc.append(self.LinLayer(self.tdim,self.net[i])) + self.tc.append(self.LinLayer(self.tdim, self.tdim)) + ''' + self.finalLayer = nn.Sequential(nn.Linear(self.tdim, self.ydim)) + + + self.X_embed = nn.Linear(self.xdim, self.tdim) + #self.T_embed = nn.Linear(self.ydim, self.tdim) + + + def LinLayer(self, dimi, dimo): + + return nn.Sequential(nn.Linear(dimi,dimo), + nn.SiLU(), + nn.LayerNorm(dimo), + nn.Dropout(p=0.1)) + + def forward(self, x): + x = self.X_embed(x) + + res_x = x + + for i in range(0,len(self.fc)): + x = self.fc[i](x) + + x = torch.add(x,res_x) + x = self.finalLayer(x) + + return x + + + + + +''' +================================================================= +Diffusion Functions +================================================================= +''' + +def timestep_embedding(timesteps, dim, max_period=10000, device=torch.device('cpu')): + """ + From https://github.com/rotot0/tab-ddpm + + Create sinusoidal timestep embeddings. + + :param timesteps: a 1-D Tensor of N indices, one per batch element. + These may be fractional. + :param dim: the dimension of the output. + :param max_period: controls the minimum frequency of the embeddings. + :return: an [N x dim] Tensor of positional embeddings. + """ + half = dim // 2 + freqs = torch.exp( + -math.log(max_period) * torch.arange(start=0, end=half, dtype=torch.float32) / half + ).to(device=timesteps.device) + args = timesteps[:, None].float() * freqs[None] + embedding = torch.cat([torch.cos(args), torch.sin(args)], dim=-1) + if dim % 2: + embedding = torch.cat([embedding, torch.zeros_like(embedding[:, :1])], dim=-1,device=device) + return embedding + +def generate_performance_weights(num_samples,num_metrics, gen_type='random'): + + weights = np.zeros((num_samples,num_metrics)) + + if gen_type == 'random': + for i in range(0,num_samples): + a = np.random.rand(1,num_metrics) + weights[i] = a/np.sum(a) + + elif gen_type == 'uniform': + samples = [] + + steps = np.linspace(0.0,1.0,11) + + for i in range(0, len(steps)): + for j in range(0,len(steps)-i): + samples.append([steps[i],steps[j],1.0-steps[i]-steps[j]]) + samples = np.array(samples) + + L = len(samples) + + print(L) + + A = np.random.randint(0,L,num_samples) + + for i in range(0,num_samples): + weights[i] = samples[A[i]] + + return weights + + + +# Now lets make a Denoise Model: + +class Denoise_ResNet_Model(torch.nn.Module): + def __init__(self, DDPM_Dict): + nn.Module.__init__(self) + + self.xdim = DDPM_Dict['xdim'] + self.ydim = DDPM_Dict['ydim'] + self.tdim = DDPM_Dict['tdim'] + self.cdim = DDPM_Dict['cdim'] + self.net = DDPM_Dict['net'] + + self.fc = nn.ModuleList() + + self.fc.append(self.LinLayer(self.tdim,self.net[0])) + + for i in range(1, len(self.net)): + self.fc.append(self.LinLayer(self.net[i-1],self.net[i])) + + self.fc.append(self.LinLayer(self.net[-1], self.tdim)) + + + self.finalLayer = nn.Sequential(nn.Linear(self.tdim, self.xdim)) + + + self.X_embed = nn.Linear(self.xdim, self.tdim) + + + self.Con_embed = nn.Sequential( + nn.Linear(self.cdim, self.tdim), + nn.SiLU(), + nn.Linear(self.tdim, self.tdim)) + + + + self.time_embed = nn.Sequential( + nn.Linear(self.tdim, self.tdim), + nn.SiLU(), + nn.Linear(self.tdim, self.tdim)) + + + def LinLayer(self, dimi, dimo): + + return nn.Sequential(nn.Linear(dimi,dimo), + nn.SiLU(), + nn.BatchNorm1d(dimo), + nn.Dropout(p=0.1)) + + + + def forward(self, x, cons, timesteps): + + + x = self.X_embed(x) + self.time_embed(timestep_embedding(timesteps, self.tdim)) + self.Con_embed(cons) + res_x = x + + for i in range(0,len(self.fc)): + x = self.fc[i](x) + + x = torch.add(x,res_x) + + x = self.finalLayer(x) + + return x + + + + +''' +============================================================================== +EMA - Exponential Moving Average: Helps with stable training +======================================================================== +EMA class from: https://github.com/azad-academy/denoising-diffusion-model/blob/main/ema.py + +''' +# Exponential Moving Average Class +# Orignal source: https://github.com/acids-ircam/diffusion_models + + +class EMA(object): + def __init__(self, mu=0.999): + self.mu = mu + self.shadow = {} + + def register(self, module): + for name, param in module.named_parameters(): + if param.requires_grad: + self.shadow[name] = param.data.clone() + + def update(self, module): + for name, param in module.named_parameters(): + if param.requires_grad: + self.shadow[name].data = (1. - self.mu) * param.data + self.mu * self.shadow[name].data + + def ema(self, module): + for name, param in module.named_parameters(): + if param.requires_grad: + param.data.copy_(self.shadow[name].data) + + def ema_copy(self, module): + module_copy = type(module)(module.config).to(module.config.device) + module_copy.load_state_dict(module.state_dict()) + self.ema(module_copy) + return module_copy + + def state_dict(self): + return self.shadow + + def load_state_dict(self, state_dict): + self.shadow = state_dict + + + +''' +======================================================================= +Trainer class modified from Tab-ddpm paper code with help from hugging face +===================================================================== +''' +class GuidedDiffusionEnv: + #def __init__(self, DDPM_Dict, Class_Dict, Reg_Dict, X,): + def __init__(self, DDPM_Dict, Class_Dict, Drag_Reg_Dict,LOA_wBulb_Reg_Dict, WL_Reg_Dict,Vol_Reg_Dict, X, X_neg, VolVec, BOAVec, DdVec): + self.DDPM_Dict = DDPM_Dict + self.datalength = self.DDPM_Dict['datalength'] + self.batch_size = self.DDPM_Dict['batch_size'] + self.Class_Dict = Class_Dict + self.Drag_Reg_Dict = Drag_Reg_Dict + self.LOA_wBulb_Reg_Dict = LOA_wBulb_Reg_Dict + self.WL_Reg_Dict = WL_Reg_Dict + self.Vol_Reg_Dict = Vol_Reg_Dict + + + self.device =torch.device(self.DDPM_Dict['device_name']) + + #Build the Diffusion Network + self.diffusion = Denoise_ResNet_Model(self.DDPM_Dict) + + + #Build Classifier Network + self.classifier = Classifier_Model(self.Class_Dict) + + #Build Regression Networks: + #self.load_trained_Drag_regressor() + + + #self.num_regressors = self.Reg_Dict['num_regressors'] + #self.load_trained_regressors() + + self.diffusion.to(self.device) + self.classifier.to(self.device) + + self.gamma = self.DDPM_Dict['gamma'] + self.lam = self.DDPM_Dict['lambda'] + + + ''' + for i in range(0,self.num_regressors): + self.regressors[i].to(self.device) + + self.dataLength = self.DDPM_Dict['datalength'] + self.batch_size = self.DDPM_Dict['batch_size'] + + self.lambdas = np.array(self.DDPM_Dict['lambdas']) + + ''' + self.data_norm = Data_Normalizer(np.array(self.DDPM_Dict['X_LL']),np.array(self.DDPM_Dict['X_UL']),self.datalength) + + # Set Up Design Data + self.X = self.data_norm.fit_Data(X) + self.X = torch.from_numpy(self.X.astype('float32')) + + #Set Up Negative Design Data + self.X_neg = self.data_norm.transform_Data(X_neg) + self.X_neg = torch.from_numpy(self.X_neg.astype('float32')) + + #Set Up Feasibility Labels + self.Cons = torch.from_numpy(np.zeros((len(self.X),1)).astype('float32')) + self.Cons_neg = torch.from_numpy(np.ones((len(self.X_neg),1)).astype('float32')) + + self.num_WL_Steps = len(VolVec[0]) #should be 101 + + + self.T_range = [.25,.67] + + self.T_vec = np.linspace(0,1,self.num_WL_Steps) + self.VolVec = VolVec + self.BOAVec = BOAVec + self.DdVec = DdVec + + + + ''' + self.X_neg = self.data_norm.transform_Data(X_neg) + + #X and Y are numpy arrays - convert to tensors + self.X_neg = torch.from_numpy(self.X_neg.astype('float32')) + self.Y = torch.from_numpy(Y.astype('float32')) + + self.Cons = torch.from_numpy(Cons.astype('float32')) + + self.Cons_neg = torch.from_numpy(Cons_neg.astype('float32')) + + + + self.X_neg = self.X_neg.to(self.device) + self.Y = self.Y.to(self.device) + ''' + + self.X = self.X.to(self.device) + self.X_neg = self.X_neg.to(self.device) + self.Cons = self.Cons.to(self.device) + self.Cons_neg = self.Cons_neg.to(self.device) + + + self.eps = 1e-8 + + self.ema = EMA(0.99) + self.ema.register(self.diffusion) + + + #set up optimizer + self.timesteps = self.DDPM_Dict['Diffusion_Timesteps'] + self.num_diffusion_epochs = self.DDPM_Dict['Training_Epochs'] + + #self.num_classifier_epochs = self.Class_Dict['Training_Epochs'] + #self.num_regressor_epochs = self.Reg_Dict['Training_Epochs'] + + lr = self.DDPM_Dict['lr'] + self.init_lr = lr + weight_decay = self.DDPM_Dict['weight_decay'] + + self.optimizer_diffusion = torch.optim.AdamW(self.diffusion.parameters(), lr=lr, weight_decay=weight_decay) + self.optimizer_classifier = torch.optim.AdamW(self.classifier.parameters(),lr=.001, weight_decay=weight_decay) + #self.optimizer_regressors = [torch.optim.AdamW(self.regressors[i].parameters(),lr=.001, weight_decay=weight_decay) for i in range(0,self.Reg_Dict['num_regressors'])] + + + + self.log_every = 100 + self.print_every = 5000 + self.loss_history = [] + + + + #Set up alpha terms + self.betas = torch.linspace(0.001, 0.2, self.timesteps).to(self.device) + + #self.betas = betas_for_alpha_bar(self.timesteps, lambda t: np.cos((t + 0.008) / 1.008 * np.pi / 2) ** 2,) + #self.betas = torch.from_numpy(self.betas.astype('float32')).to(self.device) + + self.alphas = 1. - self.betas + + self.log_alpha = torch.log(self.alphas) + self.log_cumprod_alpha = np.cumsum(self.log_alpha.cpu().numpy()) + + self.log_cumprod_alpha = torch.tensor(self.log_cumprod_alpha,device=self.device) + + + self.alphas_cumprod = torch.cumprod(self.alphas, axis=0) + self.alphas_cumprod_prev = F.pad(self.alphas_cumprod[:-1],[1,0],'constant', 0) + + self.sqrt_alphas_cumprod = torch.sqrt(self.alphas_cumprod) + self.sqrt_one_minus_alphas_cumprod = torch.sqrt(1.0 - self.alphas_cumprod) + self.sqrt_recip_alphas_cumprod = torch.sqrt(1.0 / self.alphas_cumprod) + self.sqrt_recipm1_alphas_cumprod = torch.sqrt(1.0 / self.alphas_cumprod - 1) + + self.posterior_variance = self.betas * (1. - self.alphas_cumprod_prev) / (1. - self.alphas_cumprod) + + self.sqrt_recip_alphas = torch.sqrt(1.0 / self.alphas) + + a = torch.clone(self.posterior_variance) + a[0] = a[1] + self.posterior_log_variance_clipped = torch.log(a) + self.posterior_mean_coef1 = (self.betas * torch.sqrt(self.alphas_cumprod_prev) / (1.0 - self.alphas_cumprod)) + self.posterior_mean_coef2 = ((1.0 - self.alphas_cumprod_prev)* torch.sqrt(self.alphas) / (1.0 - self.alphas_cumprod)) + + """++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + Start the training model functions + """ + def extract(self,a, t, x_shape): + b, *_ = t.shape + t = t.to(a.device) + out = a.gather(-1, t) + while len(out.shape) < len(x_shape): + out = out[..., None] + return out.expand(x_shape) + + def _anneal_lr(self, epoch_step): + #Update the learning rate + frac_done = epoch_step / self.num_diffusion_epochs + lr = self.init_lr * (1 - frac_done) + for param_group in self.optimizer_diffusion.param_groups: + param_group["lr"] = lr + + def interp(self,A,Z,z): + # This function interpolates data to approximate A(z) given values of A(Z) + + idx = np.where(Z < z)[0][-1] + + frac = (z - Z[idx])/(Z[idx+1] - Z[idx]) + + return A[idx] + frac*(A[idx+1] - A[idx]) + + def batch_train(self, batch_size=None): + ''' + This function takes in a batch of design vectors and outputs the corresponding CT values + ''' + if batch_size == None: + batch_size = self.batch_size + + A = A = np.random.randint(0,self.datalength,batch_size) + #Random Waterline + + t = np.random.uniform(self.T_range[0], self.T_range[1], (batch_size,)) + t_tens = torch.tensor(t[:,np.newaxis]).float().to(self.device) + + + #interp volume for conditioning + Vol = np.array([self.interp(self.VolVec[A[i]],self.T_vec,t[i]) for i in range(0,len(t))]) #Non-dimensionalized Waterline Length + Dd = self.DdVec[A] + BOA = self.BOAVec[A] + + x_batch = self.X[A] + + cond_batch = np.concatenate((t[:,np.newaxis],BOA[:,np.newaxis],Dd[:,np.newaxis],Vol[:,np.newaxis]),axis=1) + + cond_batch = torch.from_numpy(cond_batch.astype('float32')).to(self.device) + + return x_batch, cond_batch + + + ''' + ========================================================================= + Vanilla Diffusion + ========================================================================== + ''' + def q_sample(self,x_start, t, noise=None): + """ + qsample from https://huggingface.co/blog/annotated-diffusion + """ + if noise is None: + noise = torch.randn_like(x_start).to(self.device) + + sqrt_alphas_cumprod_t = self.extract(self.sqrt_alphas_cumprod, t, x_start.shape) + sqrt_one_minus_alphas_cumprod_t = self.extract( + self.sqrt_one_minus_alphas_cumprod, t, x_start.shape + ) + + return sqrt_alphas_cumprod_t * x_start + sqrt_one_minus_alphas_cumprod_t * noise + + + def p_loss(self,x_start, cond, t, noise=None,loss_type='l2'): + ''' + from https://huggingface.co/blog/annotated-diffusion + ''' + if noise is None: + noise = torch.randn_like(x_start) + + x_noisy = self.q_sample(x_start=x_start, t=t, noise=noise) + predicted_noise = self.diffusion(x_noisy, cond, t) + + #predicted_noise = predicted_noise.clamp(-3,3) + + if loss_type == 'l1': + loss1 = F.l1_loss(noise, predicted_noise) + elif loss_type == 'l2': + loss1 = F.mse_loss(noise, predicted_noise) + elif loss_type == "huber": + loss1 = F.smooth_l1_loss(noise, predicted_noise) + else: + raise NotImplementedError() + + return loss1 + + + + + ''' + ============================================================================== + Diffusion Training and Sampling Functions + ============================================================================== + ''' + + + def run_diffusion_step(self, x,cond): + self.optimizer_diffusion.zero_grad() + + t = torch.randint(0,self.timesteps,(self.batch_size,),device=self.device) + loss1 = self.p_loss(x,cond, t,loss_type='l2') + + loss = loss1 + loss.backward() + self.optimizer_diffusion.step() + + return loss + + def run_train_diffusion_loop(self, batches_per_epoch=100): + print('Denoising Model Training...') + self.diffusion.train() + + num_batches = self.datalength // self.batch_size + + batches_per_epoch = min(num_batches,batches_per_epoch) + + + + for i in tqdm(range(self.num_diffusion_epochs)): + + #IDX = permute_idx(self.dataLength) # get randomized list of idx for batching + + for j in range(0,batches_per_epoch): + + x_batch, cond_batch = self.batch_train() + + loss = self.run_diffusion_step(x_batch, cond_batch) + ''' + + Gaussian Diffusion (oooohhhh ahhhhhh) from TabDDPM: + ''' + #loss = self.train_step(x_batch[j]) + + self._anneal_lr(i) + + if (i + 1) % self.log_every == 0: + self.loss_history.append([i+1,float(loss.to('cpu').detach().numpy())]) + + if (i + 1) % self.print_every == 0: + print(f'Step {(i + 1)}/{self.num_diffusion_epochs} Loss: {loss}') + + + self.ema.update(self.diffusion) + #Make Loss History an np array + self.loss_history = np.array(self.loss_history) + print('Denoising Model Training Complete!') + + def fease_fn(self, x): + #From OpenAI: https://github.com/openai/guided-diffusion/blob/main/scripts/classifier_sample.py + + with torch.enable_grad(): + x_in = x.detach().requires_grad_(True) + + pred_cons = self.classifier(x_in) + + + + #log_p = torch.log(pred_cons) + + #sign = torch.sign(cons-0.5) + + grad = torch.autograd.grad(pred_cons.sum(), x_in)[0] + + #print(grad[0]) + return -grad + + def Vol_fn(self, x, cond): + #From OpenAI: https://github.com/openai/guided-diffusion/blob/main/scripts/classifier_sample.py + x_in = torch.cat((x,cond[:,0:1]),dim=1) + with torch.enable_grad(): + x_in = x_in.detach().requires_grad_(True) + + + + pred_vol = self.Vol_Reg(x_in) + + + + #log_p = torch.log(pred_cons) + + #sign = torch.sign(cons-0.5) + + grad = -2.0*(cond[:,3:4] - pred_vol)*torch.autograd.grad(pred_vol.sum(), x_in)[0] + + #print(grad[0]) + return grad[:,0:len(x[0])] + + def Drag_fn(self, x, drag_cond, g=9.81): + #From OpenAI: https://github.com/openai/guided-diffusion/blob/main/scripts/classifier_sample.py + # drag_cond = [ToD, Fn, LOA] + #concatenate the drag conditioning to the design vector + + #LOA = drag_cond[:,2:3]/ self.LOA_wBulb_Reg(x) #Calculate the desired LOA considering the bulbs + LOA = drag_cond[:,2:3] + + x_in = torch.cat((x,drag_cond[:,0:1]),dim=1) #concatenate ToD for WL Prediction + + Fn_cond = drag_cond[:,1:2]/torch.sqrt(g* LOA*self.WL_Reg(x_in)) #Calculate Froude Number for embedding + + #print(x.shape) + #print(Fn_cond.shape) + #print(LOA.shape) + + x_in = torch.cat((x_in,Fn_cond, torch.log10(LOA)),dim=1) #concatenate for drag prediction + + + + with torch.enable_grad(): + x_in = x_in.detach().requires_grad_(True) + + + + perf = self.Drag_Reg(x_in) + + + grad = torch.autograd.grad(perf.sum(), x_in)[0] + + #print(grad[0]) + return grad[:,0:len(x[0])] + + @torch.no_grad() + def p_sample(self, x, t, cons): + + time= torch.full((x.size(dim=0),),t,dtype=torch.int64,device=self.device) + + X_diff = self.diffusion(x, cons, time) + + + betas_t = self.extract(self.betas, time, x.shape) + + sqrt_one_minus_alphas_cumprod_t = self.extract( + self.sqrt_one_minus_alphas_cumprod, time, x.shape + ) + + sqrt_recip_alphas_t = self.extract(self.sqrt_recip_alphas, time, x.shape) + + """ + Compute the mean for the previous step, given a function cond_fn that + computes the gradient of a conditional log probability with respect to + x. In particular, cond_fn computes grad(log(p(y|x))), and we want to + condition on y. + + This uses the conditioning strategy from Sohl-Dickstein et al. (2015). + """ + + # Use our model (noise predictor) to predict the mean + model_mean = sqrt_recip_alphas_t * ( + x - betas_t * X_diff/ sqrt_one_minus_alphas_cumprod_t + ) + + + posterior_variance_t = self.extract(self.posterior_variance, time, x.shape) + + + fease_grad = self.fease_fn(x) + #print(gradient.detach().to('cpu')[0]) + + if t == 0: + return model_mean + else: + + noise = torch.randn_like(x,device=self.device) + # Dot product gradient to noise + return model_mean + torch.sqrt(posterior_variance_t) * (noise*(1.0-self.gamma) + self.gamma*fease_grad.float()) + #return model_mean + torch.sqrt(posterior_variance_t) * (noise + self.gamma*fease_grad.float()) + + @torch.no_grad() + def drag_p_sample(self, x, t, geom_cons, drag_cons): + + time= torch.full((x.size(dim=0),),t,dtype=torch.int64,device=self.device) + + X_diff = self.diffusion(x, geom_cons, time) + + + betas_t = self.extract(self.betas, time, x.shape) + + sqrt_one_minus_alphas_cumprod_t = self.extract( + self.sqrt_one_minus_alphas_cumprod, time, x.shape + ) + + sqrt_recip_alphas_t = self.extract(self.sqrt_recip_alphas, time, x.shape) + + """ + Compute the mean for the previous step, given a function cond_fn that + computes the gradient of a conditional log probability with respect to + x. In particular, cond_fn computes grad(log(p(y|x))), and we want to + condition on y. + + This uses the conditioning strategy from Sohl-Dickstein et al. (2015). + """ + + # Use our model (noise predictor) to predict the mean + model_mean = sqrt_recip_alphas_t * ( + x - betas_t * X_diff/ sqrt_one_minus_alphas_cumprod_t + ) + + + posterior_variance_t = self.extract(self.posterior_variance, time, x.shape) + + + fease_grad = self.fease_fn(x) + + drag_grad = self.Drag_fn(x,drag_cons) + #print(gradient.detach().to('cpu')[0]) + + if t == 0: + return model_mean + else: + + noise = torch.randn_like(x,device=self.device) + # Dot product gradient to noise + return model_mean + torch.sqrt(posterior_variance_t) * (noise*(1.0-self.gamma) + self.gamma*fease_grad.float() - self.lam[0]*drag_grad.float()) + #return model_mean + torch.sqrt(posterior_variance_t) * (noise + self.gamma*fease_grad.float() - self.lam*drag_grad.float()) + + @torch.no_grad() + def vol_drag_p_sample(self, x, t, geom_cons, drag_cons): + + time= torch.full((x.size(dim=0),),t,dtype=torch.int64,device=self.device) + + X_diff = self.diffusion(x, geom_cons, time) + + + betas_t = self.extract(self.betas, time, x.shape) + + sqrt_one_minus_alphas_cumprod_t = self.extract( + self.sqrt_one_minus_alphas_cumprod, time, x.shape + ) + + sqrt_recip_alphas_t = self.extract(self.sqrt_recip_alphas, time, x.shape) + + """ + Compute the mean for the previous step, given a function cond_fn that + computes the gradient of a conditional log probability with respect to + x. In particular, cond_fn computes grad(log(p(y|x))), and we want to + condition on y. + + This uses the conditioning strategy from Sohl-Dickstein et al. (2015). + """ + + # Use our model (noise predictor) to predict the mean + model_mean = sqrt_recip_alphas_t * ( + x - betas_t * X_diff/ sqrt_one_minus_alphas_cumprod_t + ) + + + posterior_variance_t = self.extract(self.posterior_variance, time, x.shape) + + + fease_grad = self.fease_fn(x) + + drag_grad = self.Drag_fn(x,drag_cons) + + vol_grad = self.Vol_fn(x,geom_cons) + #print(gradient.detach().to('cpu')[0]) + + if t == 0: + return model_mean + else: + + noise = torch.randn_like(x,device=self.device) + # Dot product gradient to noise + #return model_mean + torch.sqrt(posterior_variance_t) * (noise*(1.0-self.gamma) + self.gamma*fease_grad.float() - self.lam[0]*drag_grad.float() - self.lam[1]*vol_grad.float()) + return model_mean + torch.sqrt(posterior_variance_t) * (noise*(1.0-self.gamma) + self.gamma*fease_grad.float()) - self.lam[0]*drag_grad.float() - self.lam[1]*vol_grad.float() + + + @torch.no_grad() + def gen_cond_samples(self, cons): + #COND is a numpy array of the conditioning it is shape (num_samples,conditioning terms) + num_samples = len(cons) + + cons = torch.from_numpy(cons.astype('float32')) + cons = cons.to(self.device) + + #print(num_samples) #should be 1 + + x_gen = torch.randn((num_samples,self.diffusion.xdim),device=self.device) + + self.diffusion.eval() + self.classifier.eval() + + + for i in tqdm(range(self.timesteps - 1, 0, -1)): + + + x_gen = self.p_sample(x_gen, i,cons) + + + output = x_gen.cpu().detach().numpy() + + + output_scaled = self.data_norm.scale_X(output) + + return output_scaled, output + + def gen_low_drag_samples(self, Geom_COND, Drag_COND): + #Geom_COND is a numpy array of the geometric conditioning. Each row is [ToD, BoL, DoL, Vol] + #Drag_COND is a numpy array of the drag conditioning. Each Row is [ToD, U [m/s], LOA [m]] + num_samples = len(Geom_COND) + + geom_cons = torch.from_numpy(Geom_COND.astype('float32')) + geom_cons = geom_cons.to(self.device) + + drag_cons = torch.from_numpy(Drag_COND.astype('float32')) + drag_cons = drag_cons.to(self.device) + + #print(num_samples) #should be 1 + + x_gen = torch.randn((num_samples,self.diffusion.xdim),device=self.device) + + self.diffusion.eval() + self.classifier.eval() + + + + guidance_step = 16 + + for i in tqdm(range(self.timesteps - 1, guidance_step, -1)): + + + x_gen = self.drag_p_sample(x_gen, i,geom_cons, drag_cons) + + for i in tqdm(range(guidance_step, 0, -1)): + + + x_gen = self.p_sample(x_gen, i,geom_cons) + + + output = x_gen.cpu().detach().numpy() + + + output_scaled = self.data_norm.scale_X(output) + + #LOA = drag_cons[:,2:3]/ self.LOA_wBulb_Reg(x_gen) #Calculate the desired LOA considering the bulbs + OA = drag_cons[:,2:3] + LOA = LOA.cpu().detach().numpy() + + output_scaled = np.concatenate((LOA,output_scaled),axis=1) + + return output_scaled, output + + def gen_vol_drag_guided_samples(self, Geom_COND, Drag_COND): + #Geom_COND is a numpy array of the geometric conditioning. Each row is [ToD, BoL, DoL, Vol] + #Drag_COND is a numpy array of the drag conditioning. Each Row is [ToD, U [m/s], LOA [m]] + num_samples = len(Geom_COND) + + geom_cons = torch.from_numpy(Geom_COND.astype('float32')) + geom_cons = geom_cons.to(self.device) + + drag_cons = torch.from_numpy(Drag_COND.astype('float32')) + drag_cons = drag_cons.to(self.device) + + #print(num_samples) #should be 1 + + x_gen = torch.randn((num_samples,self.diffusion.xdim),device=self.device) + + self.diffusion.eval() + self.classifier.eval() + + guidance_step = 32 + + for i in tqdm(range(self.timesteps - 1, guidance_step, -1)): + + + x_gen = self.vol_drag_p_sample(x_gen, i,geom_cons, drag_cons) + + for i in tqdm(range(guidance_step, 0, -1)): + + + x_gen = self.p_sample(x_gen, i,geom_cons) + + + output = x_gen.cpu().detach().numpy() + + + output_scaled = self.data_norm.scale_X(output) + + #LOA = drag_cons[:,2:3]/ self.LOA_wBulb_Reg(x_gen) #Calculate the desired LOA considering the bulbs + LOA = drag_cons[:,2:3] + + LOA = LOA.cpu().detach().numpy() + + output_scaled = np.concatenate((LOA,output_scaled),axis=1) + + return output_scaled, output + + + ''' + ============================================================================== + Classifier and Regression Training Functions + ============================================================================== + ''' + + + def run_classifier_step(self,x,cons): + + self.optimizer_classifier.zero_grad() + + + predicted_cons = self.classifier(x) + + loss = F.binary_cross_entropy(predicted_cons, cons) #F.mse_loss(predicted_cons, cons) #F.binary_cross_entropy(predicted_cons, cons) + loss.backward() + self.optimizer_classifier.step() + + return loss + + def run_train_classifier_loop(self, batches_per_epoch=100): + + X = torch.cat((self.X,self.X_neg)) + C = torch.cat((self.Cons,self.Cons_neg)) + + test = np.random.randint(0,len(X),int(len(X)*.25)) + + X_train = X[~test] + C_train = C[~test] + + X_test = X[test] + C_test = C[test] + + print(C_train.shape) + + datalength = X_train.shape[0] + + print('Classifier Model Training...') + self.classifier.train() + + num_batches = datalength // self.batch_size + + batches_per_epoch = min(num_batches,batches_per_epoch) + + + for i in tqdm(range(self.classifier.epochs)): + + #IDX = permute_idx(self.dataLength) # get randomized list of idx for batching + + for j in range(0,batches_per_epoch): + + A = np.random.randint(0,datalength,self.batch_size) + x_batch = X_train[A] + #y_batch[j] = self.Y[IDX[j*self.batch_size:(j+1)*self.batch_size]] + cons_batch = C_train[A] + #cons_batch[j] = self.Cons[IDX[j*self.batch_size:(j+1)*self.batch_size]] + + loss = self.run_classifier_step(x_batch,cons_batch) + ''' + for i in tqdm(range(0,self.num_classifier_epochs)): + loss = self.run_classifier_step(X,C) + ''' + self.classifier.eval() + + + + + C_pred = self.classifier(X_test) + + + C_pred = C_pred.to(torch.device('cpu')).detach().numpy() + + #print(C_pred.shape) + C_pred = np.rint(C_pred) #Make it an iteger guess + + C_test = C_test.to(torch.device('cpu')).detach().numpy() + + F1 = f1_score(C_test,C_pred) + + print('F1 score: ' + str(F1)) + + print('Classifier Training Complete!') + + def Predict_Drag(self,x,drag_cond, rho=1025.0, g=9.81): + #x is a numpy array of the normalized design vector (alread) + #drag_cond is a numpy array of the drag conditioning. Each Row is [ToD, U[m/s], LOA] + + x = torch.from_numpy(x.astype('float32')).to(self.device) + + drag_cond = torch.from_numpy(drag_cond.astype('float32')).to(self.device) + + #LOA = drag_cond[:,2:3]/ self.LOA_wBulb_Reg(x) #Calculate the desired LOA considering the bulbs + + LOA = drag_cond[:,2:3] + + x = torch.cat((x,drag_cond[:,0:1]),dim=1) #concatenate ToD for WL Prediction + + Fn_cond = drag_cond[:,1:2]/torch.sqrt(g* LOA*self.WL_Reg(x)) #Calculate Froude Number for embedding + + x = torch.cat((x,Fn_cond, torch.log10(LOA)),dim=1) #concatenate for drag prediction + + self.Drag_Reg.eval() + + CT = 10**self.Drag_Reg(x) + + Drag = CT*0.5*rho*(drag_cond[:,1:2]**2)*LOA**2 #Calculate Drag: RT = CT*0.5*rho*U^2*LOA^2 + + Drag = Drag.to(torch.device('cpu')).detach().numpy() + + return Drag + + + ''' + ============================================================================== + Saving and Loading Model Functions + ============================================================================== + ''' + + def load_trained_diffusion_model(self,PATH): + #PATH is full path to the state dictionary, including the file name and extension + self.diffusion.load_state_dict(torch.load(PATH, map_location=self.device)) + + def Load_Dict(PATH): + #returns the dictionary for the DDPM_Dictionary to rebuild the model + #PATH is the path including file name and extension of the json file that stores it. + f = open(PATH) + return json.loads(f) + + + def Save_diffusion_model(self,PATH,name): + ''' + PATH is the path to the folder to store this in, including '/' at the end + name is the name of the model to save without an extension + ''' + torch.save(self.diffusion.state_dict(), PATH+name+'_diffusion.pth') + + JSON = json.dumps(self.DDPM_Dict) + f = open(PATH+name+'.json', 'w') + f.write(JSON) + f.close() + + def load_trained_classifier_model(self,PATH): + #PATH is full path to the state dictionary, including the file name and extension + self.classifier.load_state_dict(torch.load(PATH, map_location=self.device)) + + def load_trained_Drag_regression_models(self,PATH): + #label = self.Drag_Reg_Dict['Model_Paths'] + + self.Drag_Reg = Drag_Regression_ResNet(self.Drag_Reg_Dict) + self.Drag_Reg.load_state_dict(torch.load(PATH[0],map_location=self.device)) + self.Drag_Reg.to(self.device) + + self.LOA_wBulb_Reg = Regression_ResNet(self.LOA_wBulb_Reg_Dict) + self.LOA_wBulb_Reg.load_state_dict(torch.load(PATH[1],map_location=self.device)) + self.LOA_wBulb_Reg.to(self.device) + + self.WL_Reg = Regression_ResNet(self.WL_Reg_Dict) + self.WL_Reg.load_state_dict(torch.load(PATH[2],map_location=self.device)) + self.WL_Reg.to(self.device) + + self.Vol_Reg = Regression_ResNet(self.Vol_Reg_Dict) + self.Vol_Reg.load_state_dict(torch.load(PATH[3],map_location=self.device)) + self.Vol_Reg.to(self.device) + self.Drag_Reg.eval() + self.LOA_wBulb_Reg.eval() + self.WL_Reg.eval() + self.Vol_Reg.eval() + + + + + + + def Save_classifier_model(self,PATH,name): + ''' + PATH is the path to the folder to store this in, including '/' at the end + name is the name of the model to save without an extension + ''' + torch.save(self.classifier.state_dict(), PATH+name+ '.pth') + + JSON = json.dumps(self.Class_Dict) + f = open(PATH+name+ '.json', 'w') + f.write(JSON) + f.close() + + def Save_regression_models(self,PATH): + ''' + PATH is the path to the folder to store this in, including '/' at the end + + ''' + for i in range(0,len(self.regressors)): + torch.save(self.regressors[i].state_dict(), PATH + self.Reg_Dict['Model_Labels'][i] +'.pth') + + JSON = json.dumps(self.Reg_Dict) + f = open(PATH + '_regressor_Dict.json', 'w') + f.write(JSON) + f.close() diff --git a/ShipGen/tools/HullParameterization.py b/ShipGen/tools/HullParameterization.py new file mode 100644 index 0000000..fb74fee --- /dev/null +++ b/ShipGen/tools/HullParameterization.py @@ -0,0 +1,2230 @@ + #!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Mar 29 11:50:30 2022 + +@author: shannon + +Code Written by Noah Bagazinski +============================================================================================= + +The Target of this Script is to Define the functions and parameters to define a ship hull. + +This Code is a continuation of Noah Bagazinski's work in ship hull genneration using a +parametric design scheme. + +============================================================================================= +The hull parameterization is defined in five chunks: + +1) General Hull Form. + + This includes metrics such as the beam at the deck, the + hull taper, and the stern taper. + +2) The parallel mid body cross section. + + This includes parameters for the flare, chine radius, deadrise, and keel radius + +3) The Bow Form. + + This includes functions that define the bow rise (rake), the profile drift angle with respect to depth, + the profile of rocker at the front of the ship, and the profile of the location where the full breadth + is achieved for a given depth. + +4) The Stern Form + + This includes functions that define the stern profile, the convergence angle (slope of the ship at the transom), the transom cross section, the profile + of the rocker at the stern of the ship, the profile where the full beadth of the ship is achieved, and the curvature of the hull along the stern. + +5) Bulb Forms + + Bulb forms slightly modified from the parameterization defined by: + Chrismianto, D. and Kim, D.J., 2014. 'Parametric bulbous bow design using the cubic + Bezier curve and curve-plane intersection method for the minimization of ship resistance + in CFD'. Journal of Marine Science and Technology, 19(4), pp.479-492. + + functions include generation of NURBS curves to generate the meshes of the bulbous bow and stern + +""" + + +#import all the goodies: +import numpy as np +# scipy.optimize import fsolve +from matplotlib import pyplot as plt + +from stl import mesh + + + +class Hull_Parameterization: + + #Define parameters of targethull + def __init__(self, inputs): + ''' + inputs is a numpy vector that represents the parameterization of a hull. + the instantiation function generates all of the constants and factors used + in defining the parameterization of the hull. + Most of the inputs are scaled to LOA in some form + ''' + + self.LOA = inputs[0] + self.Lb = inputs[1] *self.LOA + self.Ls = inputs[2] *self.LOA + self.Bd = inputs[3]/2.0 *self.LOA#half breadth + self.Dd = inputs[4] *self.LOA + self.Bs = inputs[5] *self.Bd #half breadth, fraction of Bd + self.WL = inputs[6] *self.Dd + self.Bc = inputs[7]/2.0 *self.LOA #half breadth + self.Beta = inputs[8] + self.Rc = inputs[9]*self.Bc + self.Rk = inputs[10]*self.Dd + self.BOW = np.zeros((3,)) + self.BOW[0] = inputs[11]*0.5*self.Lb/self.Dd**2.0 + self.BOW[1] = inputs[12]*0.5*self.Lb/self.Dd + self.BK = np.zeros((2,)) + self.BK[1] = inputs[13] *self.Dd #BK_z is an input - BK_x is solved for + self.Kappa_BOW = inputs[14] + self.DELTA_BOW = np.zeros((3,)) + self.DELTA_BOW[0] = inputs[15]*0.5*self.Lb/self.Dd**2.0 + self.DELTA_BOW[1] = inputs[16]*0.5*self.Lb/self.Dd + self.DRIFT = np.zeros((3,)) + self.DRIFT[0] = inputs[17]*60.0/self.Dd**2.0 + self.DRIFT[1] = inputs[18]*60.0/self.Dd + self.DRIFT[2] = inputs[19] + self.bit_EP_S = inputs[20] + self.bit_EP_T = inputs[21] + self.TRANS = np.zeros((2,)) + self.TRANS[0] = inputs[22] + self.SK = np.zeros((2,)) + self.SK[1] = inputs[23] + self.Kappa_STERN = inputs[24] + self.DELTA_STERN = np.zeros((3,)) + self.DELTA_STERN[0] = inputs[25]*0.5*self.Ls/self.Dd**2.0 + self.DELTA_STERN[1] = inputs[26]*0.5*self.Ls/self.Dd + #self.RY_STERN = np.array(inputs[25:27]) + #self.RX_STERN = np.array(inputs[27:29]) + self.Beta_trans = inputs[27] + self.Bc_trans = inputs[28]/2.0 *self.LOA # half breadth + self.Rc_trans = inputs[29]*self.Bc_trans + self.Rk_trans = inputs[30]*self.Dd*(1-self.SK[1]) + #self.CONVERGE = np.array(inputs[33:36]) + self.bit_BB = inputs[31] + self.bit_SB = inputs[32] + self.Lbb = inputs[33] + self.Hbb = inputs[34] + self.Bbb = inputs[35] + self.Lbbm = inputs[36] + self.Rbb = inputs[37] + self.Kappa_SB = inputs[38] + self.Lsb = inputs[39] + self.HSBOA = inputs[40] + self.Hsb = inputs[41] + self.Bsb = inputs[42] + self.Lsbm = inputs[43] + self.Rsb = inputs[44] + + #Generate and Check the Forms of the Overall Hull + + self.GenGeneralHullform() + #C1 = print(self.GenralHullformConstraints()) + + self.GenCrossSection() + # C2 = print(self.CrossSectionConstraints()) + + self.GenBowForm() + #C3 = print(self.BowformConstraints()) + + self.GenSternForm() + #C4 = print(self.SternFormConstraints()) + + self.GenBulbForms() + #C5 = print(self.BulbFormConstraints()) + + + + ''' + ======================================================================= + Section 1: General Hull Form + ======================================================================= + + The General hull form is characterized by 5 characteristics: + + 0) LOA -> length overall of the vessel in [m] or = 1 + 1) Lb -> length of the bow taper in [m] or fraction of LOA + 2) Ls -> length of the stern taper in [m] or fraction of LOA + 3) Bd -> Beam at the top deck of the vessel in [m] or fraction of LOA + 4) Dd -> Depth of the vessel at the deck in [m] or fraction of LOA + 5) Bs -> Beam at the stern in [m] or fraction of LOA + 6) WL -> Waterline depts in [m] or fraction of LOA + + Constraints / NOTES to ensure realistic sizing/ shape of a hull: + 0) The length of the parallel mid body is equal to LOA-Lb-Ls = Lm + 1) Lb + Ls <= LOA + 2) Bd is not necessarily the maximum beam of the vessel. It is only the breadth of the + main deck. BOA is calculated in the Section 2: Cross Section + 3) 0 <= Bs <= BOA + 4) Lb is used to define the limits of the bow taper from the forwardmost point on the + bow rake to the point where the parallel mid-body starts. The profile of the ship + at different waterlines is dependent of the other parameters defined later in the + parameterization. + 5) Ls is used to define the limits of the stern taper from the aftmost point on the + stern rake to the point where the parallel mid-body ends. The profile of the ship + at different waterlines is dependent of the other parameters defined later in the + parameterization. + 6) WL < Dd + 7) All variables are positive or 0 + ''' + def GenGeneralHullform(self): + ''' + This funciton computes the other form factors of the general hullform + that can be calculate from the inputs + ''' + self.Lm = self.LOA - self.Ls - self.Lb + + def GenralHullformConstraints(self): + ''' + This function checks that constraints are satisfied for the hullfrom. + If no constraint violations are found, + ''' + C = np.array([-self.LOA + self.Ls+self.Lb, + self.WL - self.Dd]) + return C + ''' + ======================================================================= + Section 2: Cross Section + ======================================================================= + + The Cross Section is defined by the following inputs: + 0) Bd -> The Beam at the Deck in [m] or fraction of LOA + 1) Dd -> The Depth of the Deck in [m] or fraction of LOA + 2) Bc -> The Beam at the Chine (intersection) in [m] or fraction of LOA + 3) Dc -> The Depth of the Chine (intersection) in [m] or fraction of LOA + 4) Beta -> The deadrise angle in degrees + 5) Rc -> The Chine Radius in [m] or fraction of LOA + 6) Rk -> The keel Radius in [m] or fraction of LOA + + Constraints/ NOTES to ensure realistic sizing/ shape of a hull: + 0) 0 <= Dc < Dd + 1) Rc and Rk are agebraically limited to ensure that the radius can exist with the + given Bd,Dd,BcdC, and Beta values. + + ''' + def GenCrossSection(self): + ''' + This function calculates the constants and other form factors that will allow future + analysis of the cross section. + + ''' + + #(y,z) pair for center of keel radius + self.Rk_Center = np.array([-self.Rk*(0.5 - 0.5*np.sign(self.Rk)), + self.Rk*(0.5 + 0.5*np.sign(self.Rk))]) + #(y,z) pair for intersection of keel radius and LG line at the transom + self.Rk_LG_int = np.array([self.Rk_Center[0] + self.Rk*np.sin(np.pi*self.Beta/180.0), + self.Rk_Center[1] - self.Rk*np.cos(np.pi*self.Beta/180.0)]) + + + #solve for the lower gunwhale line: A*z + B*y + C = 0 + A = np.array([[1.0, 1.0, 1.0], + [self.Rk_LG_int[1], self.Rk_LG_int[0], 1.0], + [-(self.Rk_LG_int[0]-self.Rk_Center[0]), (self.Rk_LG_int[1]-self.Rk_Center[1]), 0.0]]) + b = np.array([1.0, 0.0, 0.0]) + + self.LG = np.linalg.solve(A,b) + + del A, b + + self.Dc = -(self.LG[1]*self.Bc + self.LG[2])/self.LG[0] + + # Upper Gunwhale Line: A*z + B*y + C = 0, where UG = [A,B,C] + A = np.array([[self.Dc, self.Bc, 1.0], + [self.Dd, self.Bd, 1.0], + [1.0, 1.0, 1.0]]) + + b = np.array([0.0,0.0,1.0]) + + self.UG = np.linalg.solve(A,b) + + del A, b + + # Calculate terms for the half beam of the cross section of the transom: + self.Rc_Center = np.zeros((2,)) #(y,z) pair for center of chine radius at the transom + self.Rc_UG_int = np.zeros((2,)) #(y,z) pair for intersection of chine radius and UG line at the transom + self.Rc_LG_int = np.zeros((2,)) #(y,z) pair for intersection of chine radius and LG line at the transom + + #make math more readable to solve the chine + A1 = self.UG[0] + B1 = self.UG[1] + theta = np.arctan2(-B1,A1) + + + if theta < 0.0: + theta = theta + np.pi + + beta = self.Beta*np.pi/180.0 + A2 = self.LG[0] + B2 = self.LG[1] + + + A = np.array([[B1, A1, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, B2, A2, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, -1.0, 0.0], + [0.0, -1.0, 0.0, 0.0, 0.0, 1.0], + [0.0, 0.0, 1.0, 0.0, -1.0, 0.0], + [0.0, 0.0, 0.0, -1.0, 0.0, 1.0]]) + + b = np.array([-self.UG[2], + -self.LG[2], + self.Rc*np.sin(theta), + self.Rc*np.cos(theta), + self.Rc*np.sin(beta), + self.Rc*np.cos(beta)]) + + C = np.linalg.solve(A,b) + + self.Rc_UG_int = C[0:2] + self.Rc_LG_int = C[2:4] + self.Rc_Center = C[4:6] + + + def CrossSectionConstraints(self): + + C = [-self.Rc_UG_int[1] + self.Dc, + -self.Rc, + -self.Bc, + -self.Dc, + self.Rc_LG_int[0] - self.Bc, + self.Rk_LG_int[0] - self.Rc_LG_int[0], + 0.00000001 -np.abs(self.Rk)] + return C + + def halfBeam_MidBody(self, z): + # This funtion calculates the half beam of the cross section at a given height, z + # If 0 > z or Dd < z, then the function returns -1 as an error + + if z < 0.0 or z > self.Dd: + return -1 + elif z >= 0.0 and z < self.Rk_LG_int[1]: + return np.sign(self.Rk)*np.sqrt((self.Rk**2) - (z-self.Rk_Center[1])**2) + self.Rk_Center[0] + elif z >= self.Rk_LG_int[1] and z < self.Rc_LG_int[1]: + return -(self.LG[0] * z + self.LG[2])/self.LG[1] + elif z >= self.Rc_LG_int[1] and z < self.Rc_UG_int[1]: + return np.sqrt((self.Rc**2) - (z-self.Rc_Center[1])**2) + self.Rc_Center[0] + else: + return -(self.UG[0] * z + self.UG[2])/self.UG[1] + + + def plot_MidBody_CrossSection(self): + + # Plot intersection points in blue + # Plot chine pt in green + # Plot Center of Rc and Rk in red + # half Beam(z) in black + + z = np.linspace(0.0, self.Dd, num = 200) + y = np.zeros((200,)) + for i in range(0,len(z)): + y[i] = self.halfBeam_MidBody(z[i]) + + + fig2, ax2 = plt.subplots() + ax2.axis('equal') + #plt.axis([0,10,0,10]) + + ax2.plot([self.Bd, self.Rc_UG_int[0], self.Rc_LG_int[0], self.Rk_LG_int[0], 0.0], + [self.Dd, self.Rc_UG_int[1], self.Rc_LG_int[1], self.Rk_LG_int[1], 0.0], 'o', color = 'blue') + ax2.plot([self.Rc_Center[0], self.Rk_Center[0]], [self.Rc_Center[1], self.Rk_Center[1]],'o' ,color = 'red') + ax2.plot([self.Bc], [self.Dc],'o' ,color = 'green') + ax2.plot(y,z,'-', color = 'black', linewidth = 0.75) + + + + + ''' + ======================================================================= + Section 3: Bow Form + ======================================================================= + + The Bow Form is defined by the following inputs: + 0) Dd -> The Depth of the Deck in [m] or fraction of LOA + 1) Lb -> The length of the bow taper in [m] or fraction of LOA + 2) Abow -> The z^2 term for Bow(z) that defines the profile of the bowrise + 3) Bbow -> The z term for Bow(z) that defines the profile of the bowrise + 4) BK_z -> The Z Point of the intersection of the Bow rise and keel rise as percentage of Dd + 5) Kappa_bow-> The X position where the Keel rise begins. percentage of Lb + 6) Adel -> z^2 term for delta(z), the x position where the max Beam is achieved for a given height + 7) Bdel -> z term for delta(z), the x position where the max Beam is achieved for a given height + 8) Adrft-> z^2 term for drift(z), the drift angle along the bowrise and keel rise + 9) Bdrft-> z term for drift(z), the drift angle along the bowrise and keel rise + 10) Cdrft-> const term for drift(z), the drift angle along the bowrise and keel rise + + These Parameters solve for 4 functions: + 0) Bow(z) -> gives the X position of the bow rise in the form Az^2 + Bz + C + 1) Keel_BOW(x) -> gives the z height of the keel rise with respect to X in the form A*(X-Kappa_BOW*Lb)^2 + 2) Delta_BOW(z) -> gives the x position between 0 and Lb where the full breadth is achieved for a given z: A(z/Dd)^2 + B(z/Dd) + C = x/Lb + 3) Drift(z) -> gives the drift angle of the bow for a given z: Az^2 + Bz + C + + These four functions define the following curve for each z: + halfBeam_Bow(x) = Y(x) = A*x^3 + Bx^2 + Cx + D for all z between 0 and Dd + Since we know two points and the derivatives of those two points + + Constraints/ NOTES to ensure realistic sizing/ shape of a hull: + 0) Kappa_BOW*Lb < delta(z=0) + 1) 0 < drift(z) < 90 for 0 <= z <= Dd (only need to check at z = 0, Dd, and -B/(2*A) if within range of z ) + 2) 0 <= BK_x < Kappa_BOW*Lb + 3) 0 <= BK_z < Dd + 4) delta(z) > Bow(z) and Keel(z) for 0 <= z <= Dd -> check z = 0,Dd,BK, Vert (Bow) and Vert (Delta) + + ''' + + def GenBowForm(self): + ''' + This funciton computes the other form factors of the Bowform + that can be calculated from the inputs + ''' + + if self.BOW[0] == 0: + Zv = -1.0 + else: + Zv = -self.BOW[1]/(2*self.BOW[0]) #Find Z of vertex of bowrise(z) + + C = np.array([self.BOW[0]*self.Dd**2.0 + self.BOW[1]*self.Dd, #Bow rise protrusion at Deck + self.BOW[0]*self.BK[1]**2.0 + self.BOW[1]*self.BK[1], #Bow rise protrusion at Bow-keel intersection + self.BOW[0]*Zv**2.0 + self.BOW[1]*Zv]) #Bowrise protrusio at vertex of bow rise eqn + + if (Zv >= self.BK[1]*self.Dd and Zv <= self.Dd): + self.BOW[2] = -np.amin(C) # If the vertex is between the BK intersect and the Deck, then it is included in the min search + else: + self.BOW[2] = -np.amin(C[0:2]) + + + # X Position of BK intersect + self.BK[0] = self.bowrise(self.BK[1]) + + + # Calculate the Keelrise equation: it is of the form X = sqrt(Z/A) + Kappa_Bow*Lb or Z = A(X-K*Lb)**2, where self.Keel = A + self.KEEL_BOW = self.BK[1]/((self.BK[0]-self.Kappa_BOW*self.Lb)**2.0) + + + #Calculate the C for the Delta equation, where C is the constant such that max(Delta(z)) = 0 between 0 and Dd + if self.DELTA_BOW[0] == 0: + Zv = -1.0 + else: + Zv = -self.DELTA_BOW[1]/(2*self.DELTA_BOW[0]) #Find Z of vertex of Delta(z) + + C = np.array([self.DELTA_BOW[0]*self.Dd**2.0 + self.DELTA_BOW[1]*self.Dd, #BDelta at Deck + 0.0, #As is, Delta(0) = 0 + self.DELTA_BOW[0]*Zv**2.0 + self.DELTA_BOW[1]*Zv]) #Bowrise protrusion at vertex of bow rise eqn + + if (Zv >= 0.0 and Zv <= self.Dd): + self.DELTA_BOW[2] = -np.amax(C) # If the vertex is between z = 0 and the Deck, then it is included in the search + else: + self.DELTA_BOW[2] = -np.amax(C[0:2]) + + #The following funcitons return the + + def bowrise(self, z): + #returns the x position of the bowrise for a given z for BK_z <= z <= Dd + return self.BOW[0]*z**2.0 + self.BOW[1]*z + self.BOW[2] + + def keelrise_bow(self, z): + #returns the x position of the keelrise at the bow for a given z for 0 <= z <= Bk_z + return -np.sqrt(z/self.KEEL_BOW) + self.Kappa_BOW*self.Lb + + def delta_bow(self, z): + #returns the x position where the full cross section width is achieved for a given z for 0 <= z <= Dd + return self.Lb + self.DELTA_BOW[0]*z**2.0 + self.DELTA_BOW[1]*z + self.DELTA_BOW[2] + + + def drift(self, z): + #returns the drift angle in radians + return np.pi*(self.DRIFT[0]*z**2.0 + self.DRIFT[1]*z + self.DRIFT[2])/180.0 + + def solve_waterline_bow(self,z): + #this function solves for a cubic function: y(half beam) = Ax^3 + Bx^2 + CX + D for the half beam of the profile between the bow/keel rise and delta for a given z for 0 <= z <= Dd + + X1 = self.bow_profile(z) + + X2 = self.delta_bow(z) + + Y2 = self.halfBeam_MidBody(z) + + A = np.array([[X1**3.0, X1**2.0, X1, 1.0], + [3.0*X1**2.0, 2*X1, 1.0, 0.0], + [X2**3.0, X2**2.0, X2, 1.0], + [3.0*X2**2.0, 2.0*X2, 1.0, 0.0]]) + + b = np.array([0.0, + np.tan(self.drift(z)), + Y2, + 0.0]) + return np.linalg.solve(A,b) + + def bow_profile(self, z): + # This assumes that z >= 0 and z <= Dd + + if z <= self.BK[1]: + X1 = self.keelrise_bow(z) + else: + X1 = self.bowrise(z) + return X1 + + def halfBeam_Bow(self, x, PROF): + #returns the halfbeam along the bow taper between the bow/keel rise and delta(z), PROF is the output of solve)waterline_bow(z) + #x is a vector + y = np.zeros((len(x),)) + for i in range(0,len(x)): + y[i] = PROF[0]*x[i]**3.0 + PROF[1]*x[i]**2.0 + PROF[2]*x[i] + PROF[3] + return y + + def bow_dydx(self, x, PROF): + #returns slope dydx of the bow taper at a height z that is defined by PROF + #x is a vecotr and function returns the vector of dydx + + dydx = np.zeros((len(x),)) + + for i in range(0,len(x)): + dydx[i] = 3.0*PROF[0]*x[i]**2.0 + 2.0*PROF[1]*x[i] + PROF[2] + + return dydx + + + def gen_waterline_bow(self, z, NUM_POINTS = 100, X = [0,1], bit_spaceOrGrid = 1): + ''' + This fuction generates a set of points [[X1,Y1] .... [X2,Y2]] that detail the curvature of the bow taper for a given z, for 0 <= z <= Dd + + it can either be created as with a number of set of points, or an even spacing based on the global x spacing (better for station plotting) + + BOOL_PTS_OR_SPACE controls whether a set number of points will produce the waterline (1) or the spacing vector will(0) + + ''' + + x1 = self.bow_profile(z) + + x2 = self.delta_bow(z) + + prof = self.solve_waterline_bow(z) + + #Set x based on spacing or grid + if bit_spaceOrGrid: + + x = np.linspace(x1,x2,NUM_POINTS) + XY = np.zeros((len(x),2)) + + else: + x = [i for i in X if (i > x1 and i <= x2)] + + + x = np.concatenate(([x1],x)) + XY = np.zeros((len(x),2)) + + + XY[0,:] = [x1, 0.0] + + y = self.halfBeam_Bow(x[1:], prof) + + XY[1:] = np.transpose([x[1:],y]) + + return XY + + def BowformConstraints(self): + #This fuction returns booleans if the bow constraints are satisfied as detailed above: + + #Check that the vertex (Zv) of the drift angle equation satisfies the constraint of the drift angle if + # if it lies within the bounds of 0 and Dd. If drift(z) is a line, or the vertex is outside the bounds, + # Then True is returned + if self.DRIFT[0] == 0.0: + Zv = -1.0 + else: + Zv = -self.DRIFT[1]/(2.0*self.DRIFT[0]) + + + if Zv >= 0.0 and Zv <= self.Dd: + vert_drift = [self.drift(Zv) - np.pi/2.0, + -self.drift(Zv)] + else: + vert_drift = [-1,-1] + + + #Check that Delta_Bow(z) is always greater than the leading edge of the ship (keelrise(z) and bow(z)) + #Check at z = 0, vertex of delta(z), vertex of bow(z), BKz, Dd + if self.DELTA_BOW[0] == 0.0: + Zv = -1.0 + else: + Zv = -self.DELTA_BOW[1]/ (2.0*self.DELTA_BOW[0]) + + if Zv >=0.0 and Zv <= self.Dd: + vert_delta_bow = (-self.delta_bow(Zv) + self.bow_profile(Zv)) + else: + vert_delta_bow = -1 + + + if self.BOW[0] == 0.0: + Zv = -1.0 + else: + Zv = -self.BOW[1]/ (2.0*self.BOW[0]) + + if Zv >=0.0 and Zv <= self.Dd: + vert_bow = (-self.delta_bow(Zv) + self.bow_profile(Zv)) + else: + vert_bow = -1 + + + + + C = [self.Kappa_BOW*self.Lb - self.delta_bow(0.0), + self.drift(0.0) - np.pi/2.0, + -self.drift(0.0) , + self.drift(self.Dd) - np.pi/2.0, + -self.drift(self.Dd), + vert_drift[0], + vert_drift[1], + -self.BK[0], + self.BK[0] - self.Kappa_BOW*self.Lb, + -self.BK[1], + self.BK[1] - self.Dd, + -self.delta_bow(self.Dd) + self.bow_profile(self.Dd), + -self.delta_bow(self.BK[1]) + self.BK[0], + vert_delta_bow, + vert_bow] + return C + + + ''' + ======================================================================= + Section 4: Stern Form + ======================================================================= + The Stern Form is defined by the following inputs: + 0) bit_EP_S -> defines whether the stern will be elliptical (1) or parabolic (0) below the SK intersect + 1) bit_EP_T -. Defines whether the stern will be elliptical (1) or parabolic (0) abover the SK intersect + 2) Bs -> The width of the stern at the deck of the ship in [m] or fraction of LOA + 3) Ls -> The length of the stern taper in [m] or fraction of LOA + 4) A_trans -> The A term that defines the transom slope X = Az + B + 5) SKz -> The Z Point of the intersection of the Stern rise and transom as percentage of Dd + 6) Kappa_STERN -> The X position where the Stern rise begins aft of the end of the parallel midbody as a fraction of Ls + 7) Adel -> z^2 term for delta_stern(z), the x position where the max Beam is achieved for a given height, + 8) Bdel -> z term for delta_stern(z), the x position where the max Beam is achieved for a given height + 9) Bc_trans -> The beam of the chine point at the transom in [m] or fraction of LOA + 10) Dc_trans -> The depth of the chine point at the transom in [m] or fraction of LOA + 11) Rc_trans -> The Chine radius of the chine at the transom in [m] or fraction of LOA + 12) Rk_trans -> the keel radius of the chine at the transom in [m] or fraction of LOA + + + + REMOVE THESE FOR NOW + 7) A_Ry-> z term for Ry(z), the y-raduis of the ellipse at the stern of the ship + 8) B_Ry-> const for Ry(z), the y-raduis of the ellipse at the stern of the ship + 9) A_Rx-> z term for Rx(z), the x-raduis of the ellipse at the stern of the ship + 10) B_Rx-> const for Rx(z), the x-raduis of the ellipse at the stern of the ship + + 15) AconvT -> the z^2 term for Converge Angle(z) the tangent angle of the gunwhale at the transom + 16) BconvT -> the z term for Converge Angle(z) the tangent angle of the gunwhale at the transom + 17) CconvT -> the const term for Converge Angle(z) the tangent angle of the gunwhale at the transom + + + These Parameters solve for 7 functions: + 0) Transom(z) -> gives the X position of the transom in the form Az + B + 1) Sternrise(x) -> gives the z height of the stern rise with respect to X in the form A*(X-Kappa*Ls)^2 + 2) Delta_Stern(z) -> gives the x position between LOA-Ls and LOA where the full breadth is achieved for a given z: A(z)^2 + B(z) + C = X + 3) halfBeam_transom(z) -> gives the halfbeam of the transom for z between SKz and Dd + + REMOVE THESE FOR BIW + 3) Converge(z) -> gives the convergence tangent angle of the gunwhale at the transom for a given z: Az^2 + Bz + C + 4) Ry(z) -> gives the y radius of the stern ellipse in the form Ry = Az + B + 5) Rx(z) -> gives the x radius of the stern ellipse in the form Rx = Az + B + + + These four functions define the following curve for each z: + halfBeam_Stern(x) = Y(x) = Parabola + Ellipse for all z between 0 and Dd + + Constraints/ NOTES to ensure realistic sizing/ shape of a hull: + 0) Lb+Lm + Kappa_Stern*Ls > delta_Stern(z=0) + 1) 0 < converge(z) < 90 for 0 <= z <= Dd (only need to check at z = 0, Dd, and -B/(2*A) if within range of z ) + 2) 0 <= SK_x > Lb+Lm+ Kappa*Ls + 3) 0 <= SK_z < Dd + 4) delta(z) < Transom(z) and Sternrise(z) for 0 <= z <= Dd + + + + ''' + def GenSternForm(self): + + # Recalculate SK to be a value instead of a percentage + self.SK[1] = self.SK[1]*self.Dd + + # Solve for the B value such that max(Transom(z)) = LOA + if self.TRANS[0] >= 0.0: + self.TRANS[1] = self.LOA - self.TRANS[0]*self.Dd + else: + self.TRANS[1] = self.LOA - self.TRANS[0]*self.SK[1] + + #calculate the x value for the SK intersect + self.SK[0] = self.transom(self.SK[1]) + + # find the constant term in the sternrise equation: z = A(x-Lb+Lm+Ls*Kappa_stern)^2 + self.STERNRISE = self.SK[1]/(self.SK[0] - (self.Lb + self.Lm + self.Ls*self.Kappa_STERN))**2.0 + + #Calculate the C for the Delta_stern equation, where C is the constant such that min(Delta_stern(z)) = 0 for z between 0 and Dd + if self.DELTA_STERN[0] == 0: + Zv = -1.0 + else: + Zv = -self.DELTA_STERN[1]/(2*self.DELTA_STERN[0]) #Find Z of vertex of Delta(z) + + C = np.array([self.DELTA_STERN[0]*self.Dd**2.0 + self.DELTA_STERN[1]*self.Dd, #Stern Delta at Deck + 0.0, #As is, Delta_Stern(0) = 0 + self.DELTA_STERN[0]*Zv**2.0 + self.DELTA_STERN[1]*Zv]) #vertex of Delta_STERN equation + + if (Zv >= 0.0 and Zv <= self.Dd): + self.DELTA_STERN[2] = -np.amin(C) # If the vertex is between z = 0 and the Deck, then it is included in the search + else: + self.DELTA_STERN[2] = -np.amin(C[0:2]) + + + #(y,z) pair for center of keel radius + self.Rk_Center_trans = np.array([-self.Rk_trans*(0.5 - 0.5*np.sign(self.Rk_trans)), + self.SK[1] + self.Rk_trans*(0.5 + 0.5*np.sign(self.Rk_trans))]) + #(y,z) pair for intersection of keel radius and LG line at the transom + self.Rk_LG_int_trans = np.array([self.Rk_Center_trans[0] + self.Rk_trans*np.sin(np.pi*self.Beta_trans/180.0), + self.Rk_Center_trans[1] - self.Rk_trans*np.cos(np.pi*self.Beta_trans/180.0)]) + + + #solve for the lower gunwhale line: A*z + B*y + C = 0 + A = np.array([[1.0, 1.0, 1.0], + [self.Rk_LG_int_trans[1], self.Rk_LG_int_trans[0], 1.0], + [-(self.Rk_LG_int_trans[0]-self.Rk_Center_trans[0]), (self.Rk_LG_int_trans[1]-self.Rk_Center_trans[1]), 0.0]]) + b = np.array([1.0, 0.0, 0.0]) + + self.LG_trans = np.linalg.solve(A,b) + + del A, b + + self.Dc_trans = -(self.LG_trans[1]*self.Bc_trans + self.LG_trans[2])/self.LG_trans[0] + + # Upper Gunwhale Line: A*z + B*y + C = 0, where UG = [A,B,C] + A = np.array([[self.Dc_trans, self.Bc_trans , 1.0], + [self.Dd, self.Bs, 1.0], + [1.0, 1.0, 1.0]]) + + b = np.array([0.0,0.0,1.0]) + + self.UG_trans = np.linalg.solve(A,b) + + del A, b + + # Calculate terms for the half beam of the cross section of the transom: + self.Rc_Center_trans = np.zeros((2,)) #(y,z) pair for center of chine radius at the transom + self.Rc_UG_int_trans = np.zeros((2,)) #(y,z) pair for intersection of chine radius and UG line at the transom + self.Rc_LG_int_trans = np.zeros((2,)) #(y,z) pair for intersection of chine radius and LG line at the transom + + #make math more readable to solve the chine + A1 = self.UG_trans[0] + B1 = self.UG_trans[1] + theta = np.arctan2(-B1,A1) + + if theta < 0.0: + theta = theta + np.pi + + beta = self.Beta_trans*np.pi/180.0 + A2 = self.LG_trans[0] + B2 = self.LG_trans[1] + + + A = np.array([[B1, A1, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, B2, A2, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, -1.0, 0.0], + [0.0, -1.0, 0.0, 0.0, 0.0, 1.0], + [0.0, 0.0, 1.0, 0.0, -1.0, 0.0], + [0.0, 0.0, 0.0, -1.0, 0.0, 1.0]]) + + b = np.array([-self.UG_trans[2], + -self.LG_trans[2], + self.Rc_trans*np.sin(theta), + self.Rc_trans*np.cos(theta), + self.Rc_trans*np.sin(beta), + self.Rc_trans*np.cos(beta)]) + + C = np.linalg.solve(A,b) + + self.Rc_UG_int_trans = C[0:2] + self.Rc_LG_int_trans = C[2:4] + self.Rc_Center_trans = C[4:6] + + + + + def transom(self, z): + #returns the x position of the transom for a given z fr SK_z <= z <= Dd + return self.TRANS[0]*z + self.TRANS[1] + + def sternrise(self, z): + #returns the x position of the sternrise for a given z for 0 <= z <= SK_z + return np.sqrt(z/self.STERNRISE) + self.Lb + self.Lm + self.Ls*self.Kappa_STERN + + def stern_profile(self, z): + # shows the profile of the stern without the bulbous stern: + + if self.bit_SB and z <= self.WL*self.HSBOA: + return self.SB_Prof[0] # If there is a bulbous stern, we want to form the profile to lead into the SB. + + else: + if z <= self.SK[1]: + return self.sternrise(z) + else: + return self.transom(z) + + def delta_stern(self, z): + #returns the starting position of the stern taper at a given heigt + return self.Lb + self.Lm + self.DELTA_STERN[0]* z**2.0 + self.DELTA_STERN[1]*z + self.DELTA_STERN[2] + + + def halfBeam_Transom(self, z): + #Returns the x,y pair of the transom at a height z. This assumes that SK_z <= z <= Dd. otherwise y returns -1 + x = self.stern_profile(z) + + if z <= self.SK[1] and z > 0.0: + y = 0.0 + elif z < 0.0 or z > self.Dd: + y = -1.0 + + elif z > self.SK[1] and z < self.Rk_LG_int_trans[1]: + y = np.sign(self.Rk_trans)*np.sqrt((self.Rk_trans**2) - (z-self.Rk_Center_trans[1])**2) + self.Rk_Center_trans[0] + elif z >= self.Rk_LG_int_trans[1] and z < self.Rc_LG_int_trans[1]: + y = -(self.LG_trans[0] * z + self.LG_trans[2])/self.LG_trans[1] + elif z >= self.Rc_LG_int_trans[1] and z < self.Rc_UG_int_trans[1]: + y = np.sqrt((self.Rc_trans**2) - (z-self.Rc_Center_trans[1])**2) + self.Rc_Center_trans[0] + else: + y = -(self.UG_trans[0] * z + self.UG_trans[2])/self.UG_trans[1] + + return [x,y] + + def plot_Transom_CrossSection(self): + + # Plot intersection points in blue + # Plot chine pt in green + # Plot Center of Rc and Rk in red + # half Beam(z) in black + + z = np.linspace(self.SK[1], self.Dd, num = 200) + y = np.zeros((200,2)) + for i in range(0,len(z)): + y[i] = self.halfBeam_Transom(z[i]) + + + fig1,ax1 = plt.subplots() + ax1.axis('equal') + #plt.axis([0,10,0,10]) + + ax1.plot([self.Bs, self.Rc_UG_int_trans[0], self.Rc_LG_int_trans[0], self.Rk_LG_int_trans[0], 0.0], + [self.Dd, self.Rc_UG_int_trans[1], self.Rc_LG_int_trans[1], self.Rk_LG_int_trans[1], self.SK[1]], 'o', color = 'blue') + ax1.plot([self.Rc_Center_trans[0], self.Rk_Center_trans[0]], [self.Rc_Center_trans[1], self.Rk_Center_trans[1]],'o' ,color = 'red') + ax1.plot([self.Bc_trans], [self.Dc_trans],'o' ,color = 'green') + ax1.plot(y[:,1],z,'-', color = 'black', linewidth = 0.75) + + def halfBeam_Stern(self, x, PROF): + #returns the halfbeam along the stern taper between delta(z) and stern_profile(z), PROF is the output of solve_waterline_stern(z) + # x is a vector + y = np.zeros((len(x),)) + + if PROF[0]: + for i in range(0,len(x)): + y[i] = np.sqrt( np.abs(PROF[5]**2.0 * (1.0 - ((x[i] - PROF[6])/PROF[4])**2.0))) + PROF[7] #ellipse + else: + for i in range(0,len(x)): + + y[i] = PROF[1]*x[i]**2.0 + PROF[2]*x[i] + PROF[3] #parabola + + return y + + def stern_dydx(self, x, PROF): + #returns slope dydx of the stern taper at a height z that is defined by PROF + #x is a vecotr and function returns the vector of dydx + + dydx = np.zeros((len(x),)) + + if PROF[0]: + for i in range(0,len(x)): + + dydx[i] = -PROF[5]*(x[i]-PROF[6])/(PROF[4]**2.0) * 1/np.sqrt(np.abs(1.0 - ((x[i] - PROF[6])/PROF[4])**2.0)) + + else: + for i in range(0,len(x)): + + dydx[i] = 2.0*PROF[1]*x[i] + PROF[2] + + return dydx + + + + def solve_waterline_stern(self, z): + #returns PROF, a parabola [A,B,C], an ellipse [Rx, Ry, Cx, Cy] of the two curves such that they are tangent at the intersection + # Compares bit_EP_S and bit_EP_T -> if the curve is a parabola, the ellipse values in PROF are 0 and vice versa + # PROF = [A,B,C, Rx, Ry, Cx, Cy] + + x1 = self.delta_stern(z) + y1 = self.halfBeam_MidBody(z) + + [x2,y2] = self.halfBeam_Transom(z) + + PROF = np.zeros((8,)) + + if z >= self.SK[1]: + if self.bit_EP_T: + #If the curve is at the transom and the curve is an ellipse + Rx = x2-x1 + Ry = y1-y2 + Cx = x1 + Cy = y2 + PROF[0] = 1 + PROF[4:8] = np.array([Rx,Ry,Cx,Cy]) + else: + #If the curve is at the transom and the curve is a parabola: + A = np.array([[x1**2.0, x1, 1.0], + [x2**2.0, x2, 1.0], + [2.0*x1, 1.0, 0.0]]) + + b = np.array([y1, y2, 0.0]) + + C = np.linalg.solve(A,b) + + PROF[1:4] = C + else: + if self.bit_EP_S: + #If the curve is below the transom and the curve is an ellipse + Rx = x2-x1 + Ry = y1-y2 + Cx = x1 + Cy = y2 + PROF[0] = 1 + PROF[4:8] = np.array([Rx,Ry,Cx,Cy]) + else: + #If the curve is below the transom and the curve is a parabola: + A = np.array([[x1**2.0, x1, 1.0], + [x2**2.0, x2, 1.0], + [2.0*x1, 1.0, 0.0]]) + + b = np.array([y1, y2, 0.0]) + + C = np.linalg.solve(A,b) + + PROF[1:4] = C + + return PROF + + + + def gen_waterline_stern(self, z, NUM_POINTS = 100, X = [0,1], bit_spaceOrGrid = 1): + ''' + This fuction generates a set of points [[X1,Y1] .... [X2,Y2]] that detail the curvature of the bow taper for a given z, for 0 <= z <= Dd + + it can either be created as with a number of set of points, or an even spacing based on the global x spacing (better for station plotting) + + BOOL_PTS_OR_SPACE controls whether a set number of points will produce the waterline (1) or the spacing vector will(0) + + ''' + x1 = self.delta_stern(z) + + x2 = self.stern_profile(z) + + + + prof = self.solve_waterline_stern(z) + + + + if bit_spaceOrGrid: + x = np.linspace(x1,x2,NUM_POINTS) + XY = np.zeros((len(x),2)) + + else: + x = [i for i in X if (i >= x1 and i < x2)] + x = np.concatenate((x,[x2])) + XY = np.zeros((len(x),2)) + + y = self.halfBeam_Stern(x[0:-1], prof) + + + XY[0:-1] = np.transpose([x[0:-1],y]) + + #set the last element in the array to be the transom point + XY[-1] = self.halfBeam_Transom(z) + + + return XY + + + + def SternformConstraints(self): + # this is an incomplete list of geometric constrains for the hull form + + #Check that Delta_Bow(z) is always greater than the leading edge of the ship (keelrise(z) and bow(z)) + #Check at z = 0, vertex of delta(z), vertex of bow(z), BKz, Dd + if self.DELTA_STERN[0] == 0.0: + Zv = -1.0 + else: + Zv = -self.DELTA_STERN[1]/ (2.0*self.DELTA_STERN[0]) + + if Zv >=0.0 and Zv <= self.Dd: + vert_delta_stern = (self.delta_stern(Zv) - self.stern_profile(Zv)) + else: + vert_delta_stern = -1 + + + + C = [self.delta_stern(0.0) - (self.Lb + self.Lm + self.Ls*self.Kappa_STERN), + self.delta_stern(self.SK[1]) - self.SK[0], + vert_delta_stern, + self.delta_stern(self.Dd) - self.stern_profile(self.Dd), + (self.Lb + self.Lm + self.Ls*self.Kappa_STERN) - self.SK[0], + self.Bc_trans - self.halfBeam_MidBody(self.Dc_trans), + -self.Rc_UG_int_trans[1] + self.Dc_trans, + -self.Rc_trans, + -self.Bc_trans, + -self.Dc_trans, + self.Rc_LG_int_trans[0] - self.Bc_trans, + self.Rk_LG_int_trans[0] - self.Rc_LG_int_trans[0]] + + + return C + + + + + ''' + ======================================================================= + Section 5: Bulb Forms + ======================================================================= + The Bulb Forms are defined by the following inputs: + 0) bit_BB -> Bit that defines whether there is a bublous bow (1) or not (0) + 1) bit_SB -> Bit that defines whether there is a bublous stern (1) or not (0) + 2) Lbb -> Length of the bulbous bow (BB) fwd the foward perpendicular as a fraction of LOA + 3) Hbb -> Height of widest part of the BB as a fraction of the WL + 4) Bbb -> max. Width of the BB at the FP as a fraction of Bd + 5) Lbbm -> A midpoint along the length of the BB as a fraction of Lbb -> represents the position of max beam of the BB + 6) Rbb -> radius that fillets the BB to the hull (ratio from 0 to 1) -> solved as a cubic function + 7) Kappa_SB -> Position where the Stern Bulb diverges from the hull as a fraction of Ls + 8) Lsb -> Length of the stern bulb (SB) as a fraction of LOA + 9) HSBOA -> Overall Height of the SB as a fraction of WL + 9) Hsb -> Height of widest part of the SB as a fraction of HSBOA + 10) Bsb -> max. Width of the SB at the Kappa_SB as a fraction of Bd + 11) Lsbm -> A midpoint along the length of the SB as a fraction of Lsb -> represents the position of max beam of the BB + 12) Rsb -> radius that fillets the SB to the hull (ratio from 0 to 1) -> solved as a cubic function + + + + These Parameters solve 3 functions each for the Bulbous Bow and Bulbous Stern: + 0) Outline: Definition of upper and lower ellipse that define the profile of the bulb + 1) Profile of Max width : a parabola that is tangent to an ellipse at the longitudinal mid point. + 2) Cross Section Generator: Solves for Rx and Ry of an upper and lower ellipse that solves for a cross section of the bulb + 3) + + + Constraints/ NOTES to ensure realistic sizing/ shape of a bulb: + 0) Cross Section of bulbs at Starting position need to be encompassed by Half Beam Mid Body + 1) Rk >= 0 + 2) 0.0 < (Hbb and Hsb) < 1.0 + 3) 0.0 < (Bbb and Bsb) < 1.0 + 4) 0.0 < (Lbb and Lsb) < TBD (seems a bit outrageous) (but not infeasible technically) + 5) -1.0 < (Lbbm and Lsbm) < 1.0 + 6) 0.0 < HSBOA < 1.0 + + + + + + ''' + def GenBulbForms(self): + ''' + This function generates Prof for the bulbous bow and bulbous stern. + + it is composted of the following elements in this order: + 0a) FP = forward perpendicular: the start point were the BB diverges from the hull. + 0b) SBs = start of stern bulb: the starting point where the SB diverges from the hull + 1) Rz_U = upper z radius for an ellipse that defines the top half of the bulb + 2) Rz_L = lower z radius for an ellipse that defines the bottom half of the bulb + 3) Ry -> y radius of an ellipse that defines the max width of the bulb + 4) Rx -> x Radius of leading edge of bulb + 5) Cx -> X center of aforementioned ellipse + 6) Cy -> Y center of aforementioned ellipse + + + ''' + + self.BB_Prof = np.zeros((7,)) + self.SB_Prof = np.zeros((7,)) + + if self.bit_BB: + + FP = self.bow_profile(self.WL) + + + self.BB_Prof[0] = FP + self.BB_Prof[1] = (1.0 - self.Hbb)*self.WL + self.BB_Prof[2] = self.Hbb*self.WL + self.BB_Prof[3] = self.halfBeam_MidBody(self.BB_Prof[2]) *self.Bbb + self.BB_Prof[4] = self.Lbb*self.LOA*(1.0-self.Lbbm) + self.BB_Prof[5] = FP - self.LOA*self.Lbb*self.Lbbm + self.BB_Prof[6] = 0.0 + + + if self.bit_SB: + SBs = self.Kappa_SB*self.Ls + self.Lm + self.Lb + + self.SB_Prof[0] = SBs + self.SB_Prof[1] = (1.0 - self.Hsb)*self.WL*self.HSBOA + self.SB_Prof[2] = self.Hsb*self.WL*self.HSBOA + self.SB_Prof[3] = self.halfBeam_MidBody(self.SB_Prof[2])*self.Bsb + self.SB_Prof[4] = self.Lsb*self.LOA*(1.0-self.Lsbm) + self.SB_Prof[5] = SBs + self.LOA*self.Lsb*self.Lsbm + self.SB_Prof[6] = 0.0 + + + def BB_profile(self,z): + #assumes z is between WL and 0.0 + #returns x position of leading edge of SB + if z >= self.BB_Prof[2]: + return self.BB_Prof[5] - np.sqrt(np.abs(1.0 - ((z-self.Hbb*self.WL)/self.BB_Prof[1])**2.0))*self.BB_Prof[4] + else: + return self.BB_Prof[5] - np.sqrt(np.abs(1.0 - ((z-self.Hbb*self.WL)/self.BB_Prof[2])**2.0))*self.BB_Prof[4] + + def halfBeam_BB(self, z, x): + #returns the half breadth of the BB at height z and position x. x is a vector + #assumes z is between 0 and WL assumes x greater than BB_profile(z) + + if z >= self.BB_Prof[2]: + Rz = self.BB_Prof[1] + else: + Rz = self.BB_Prof[2] + + Ry = np.sqrt(np.abs(1.0 - ((z - self.BB_Prof[2])/Rz)**2.0))*self.BB_Prof[3] + + Rx = self.BB_Prof[5] - self.BB_profile(z) + + y = np.zeros((len(x),)) + + for i in range(0,len(x)): + + if x[i] >= self.BB_Prof[5]: + y[i] = Ry + else: + y[i] = np.sqrt(np.abs(1.0 - ((x[i] - self.BB_Prof[5])/Rx)**2.0))*Ry + + return y + + + + def BB_dydx(self,z,x): + #This function computes the slope dy/dx slope of the bulbous bow at height z and position x. This assumes x is within the bulbous bow x-range + # x is a vector + + dydx = np.zeros((len(x),)) + + if z >= self.BB_Prof[2]: + Rz = self.BB_Prof[1] + else: + Rz = self.BB_Prof[2] + + Ry = np.sqrt(np.abs(1.0 - ((z - self.BB_Prof[2])/Rz)**2.0))*self.BB_Prof[3] + Rx = self.BB_Prof[5] - self.BB_profile(z) + + for i in range(0,len(x)): + + if x[i] >= self.BB_Prof[5]: + dydx[i] = 0.0 + + else: + + dydx[i] = -Ry*(x[i]-self.BB_Prof[5])/(Rx**2.0) * 1/np.sqrt(np.abs(1.0 - ((x[i] - self.BB_Prof[5])/Rx)**2.0)) + + return dydx + + + + def SB_profile(self,z): + #assumes z is between WL*HSBOA and 0.0 + #returns x position of trailing edge of SB + if z >= self.Hsb*self.WL*self.HSBOA: + # print(z) + # print((1- ((z-self.Hsb*self.WL*self.HSBOA)/self.SB_Prof[1])**2.0)) + return self.SB_Prof[5] + np.sqrt(np.abs(1.0 - ((z-self.Hsb*self.WL*self.HSBOA)/self.SB_Prof[1])**2.0))*self.SB_Prof[4] + else: + return self.SB_Prof[5] + np.sqrt(np.abs(1.0 - ((z-self.Hsb*self.WL*self.HSBOA)/self.SB_Prof[2])**2.0))*self.SB_Prof[4] + + def halfBeam_SB(self, z, x): + #returns the half breadth of the BB at height z and position x. x is a vector + #assumes z is between 0 and WL*HSBOA assumes x less than SB_profile(z) + + if z >= self.Hsb*self.WL*self.HSBOA: + Rz = self.SB_Prof[1] + else: + Rz = self.SB_Prof[2] + + Ry = np.sqrt(np.abs(1.0 - ((z - self.Hsb*self.WL*self.HSBOA)/Rz)**2.0))*self.SB_Prof[3] + + Rx = self.SB_profile(z) - self.SB_Prof[5] + + y = np.zeros((len(x),)) + + for i in range(0,len(x)): + + if x[i] <= self.SB_Prof[5]: + y[i] = Ry + else: + y[i] = np.sqrt(np.abs(1.0 - ((x[i] - self.SB_Prof[5])/Rx)**2.0))*Ry + + return y + + def SB_dydx(self,z,x): + #This function computes the slope dy/dx slope of the bulbous bow at height z and position x. This assumes x is within the bulbous bow x-range + # x is a vector + if z >= self.Hsb*self.WL*self.HSBOA: + Rz = self.SB_Prof[1] + else: + Rz = self.SB_Prof[2] + + Ry = np.sqrt(np.abs(1.0 - ((z - self.Hsb*self.WL*self.HSBOA)/Rz)**2.0))*self.SB_Prof[3] + + Rx = self.SB_profile(z) - self.SB_Prof[5] + + dydx = np.zeros((len(x),)) + + for i in range(0, len(x)): + + if x[i] <= self.SB_Prof[5]: + dydx[i] = 0.0 + + else: + dydx[i] = -Ry*(x[i]-self.SB_Prof[5])/(Rx**2.0) * 1/np.sqrt(np.abs(1.0 - ((x[i] - self.SB_Prof[5])/Rx)**2.0)) + + return dydx + + + + + + def gen_waterline_bow_BB(self, z, NUM_POINTS = 100, X = [0,1], bit_spaceOrGrid = 1): + #this function returns a set of [X,Y] points that accounts for the shape and fillet radius of a bulbous bow on the bow profile + + a = NUM_POINTS + if z >= self.WL: + #If z is above the Ship's waterline, then the bulbous bow does not exist in that section + return self.gen_waterline_bow(z, NUM_POINTS = a, X = X, bit_spaceOrGrid = bit_spaceOrGrid) + + else: + PROF = self.solve_waterline_bow(z) + + x1 = self.BB_profile(z) + + x2 = self.delta_bow(z) + + + #Set x based on spacing or grid + if bit_spaceOrGrid: + + x = np.linspace(x1,x2,NUM_POINTS) + XY = np.zeros((len(x),2)) + + else: + x = [i for i in X if (i > x1 and i <= x2)] + + x = np.concatenate(([x1],x)) + + XY = np.zeros((len(x),2)) + + + + #Find most likely point where BB intersects Bow Curve + A = PROF.copy() + + A[3] = A[3] - self.halfBeam_BB(z, [self.BB_Prof[5]]) + + + ROOTS = np.roots(A) + + x_int = np.amin(np.real([i for i in ROOTS if (i >= self.bow_profile(z) and i >= self.BB_profile(z))])) #Need to call np.real as there will be instances where 2nd and 3rd roots of PROF will be imaginary. calling np.real to clean this up + + + ''' + #ind start and ending points for Rbb -> not quite a circular radius, but it is a cubic fillet (sorta counts) + Xrad are the points forward and aft where fillet will start. + Xrad[0] is Rbb fraction of distance between interesect and fwd profie of BB at z + Xrad[1] is Rbb fraction of the distance bewtween the insect and delta_bow(z) + ''' + dx = abs(np.amin([x_int - self.BB_profile(z), self.delta_bow(z) - x_int])) #Distance over which fillet occurs # Need to add abs to avoid dumb errors + + + + Xrad = [(x_int - dx*self.Rbb/2.0), (x_int + dx*self.Rbb/2.0)] + + + # going to build quartic systems of Eqn to solve. only tricky thing: what is dydx at Xrad[0] + + Yrad = [self.halfBeam_BB(z, [Xrad[0]])[0], self.halfBeam_Bow([Xrad[1]], PROF)[0]] + dydx = [self.BB_dydx(z, [Xrad[0]])[0], self.bow_dydx([Xrad[1]], PROF)[0]] + + + ''' + Rbb is quartic systems of eqns + 5 boundary conditions: + both (Xrad, Yrad) on fillet curve + both dydx are matched at ends + dydx halfway between BCs is mean of BC dydx and avg slope of BC end points + + ''' + + #dydx_mean = (2.0*((Yrad[1] - Yrad[0])/(Xrad[1]-Xrad[0])) + dydx[0] + dydx[1])/4.0 + + Arad = np.array([[Xrad[0]**4.0, Xrad[0]**3.0, Xrad[0]**2.0, Xrad[0], 1.0], + [Xrad[1]**4.0, Xrad[1]**3.0, Xrad[1]**2.0, Xrad[1], 1.0], + [4.0*Xrad[0]**3.0, 3.0*Xrad[0]**2.0, 2.0*Xrad[0], 1.0, 0.0], + [4.0*Xrad[1]**3.0, 3.0*Xrad[1]**2.0, 2.0*Xrad[1], 1.0, 0.0], + [12.0*Xrad[0]**2.0, 6.0*Xrad[0]**2.0, 2.0, 0.0, 0.0]]) + + + + brad = np.array([Yrad[0], Yrad[1], dydx[0], dydx[1], 0.0]) #dydx_mean]) + + + PROFrad = np.linalg.solve(Arad,brad) + + XY[0,:] = [x1, 0.0] + + xbb = [i for i in x if (i > x1 and i <= Xrad[0])] + + xbbrad = [i for i in x if (i > Xrad[0] and i < Xrad[1])] + + xbow = [i for i in x if (i >= Xrad[1] and i <= x2)] + + ybb = self.halfBeam_BB(z, xbb) + + ybbrad = np.zeros((len(xbbrad),)) + + for i in range(0,len(xbbrad)): + ybbrad[i] = PROFrad[0]*xbbrad[i]**4.0 + PROFrad[1]*xbbrad[i]**3.0 + PROFrad[2]*xbbrad[i]**2.0 + PROFrad[3]*xbbrad[i] + PROFrad[4] + + ybow = self.halfBeam_Bow(xbow, PROF) + + + + XY[1:,0] = x[1:] + + XY[1:,1] = np.concatenate((ybb,ybbrad,ybow)) + + + return XY + + + def gen_waterline_stern_SB(self, z, NUM_POINTS = 100, X = [0,1], bit_spaceOrGrid = 1): + #this function returns a set of [X,Y] points that accounts for the shape and fillet radius of a bulbous bow on the bow profile + + a = NUM_POINTS + if z >= self.WL*self.HSBOA: #If z is above the Ship's waterline, then the bulbous bow does not exist in that section + return self.gen_waterline_stern(z, NUM_POINTS = a, X=X, bit_spaceOrGrid=bit_spaceOrGrid) + + else: + #Set up half beam for stern at z + PROF = self.solve_waterline_stern(z) + + x1 = self.delta_stern(z) + + x2 = self.SB_profile(z) + + #Create x distribution + if bit_spaceOrGrid: + x = np.linspace(x1,x2,NUM_POINTS) + XY = np.zeros((len(x),2)) + else: + x = [i for i in X if (i >= x1 and i < x2)] + x = np.concatenate((x,[x2])) + XY = np.zeros((len(x),2)) + + + # This y intersect is most likely place that y in + y_int = self.halfBeam_SB(z, [self.SB_Prof[5]])[0] #ad the [0] so that y_int is interpretted as a float + + if y_int >= self.halfBeam_Stern([self.delta_stern(z)], PROF): + y_int = self.halfBeam_Stern([self.delta_stern(z)], PROF) + x_int = self.SB_Prof[0] + + else: + if PROF[0]: + x_int = PROF[4]*np.sqrt(abs(1 - ((y_int-PROF[7])/PROF[5])**2.0)) + PROF[6] + else: + ROOTS = np.roots([PROF[1], PROF[2], PROF[3]-y_int]) + + x_int = np.amax(np.real([i for i in ROOTS if (i >= self.delta_stern(z))])) + + + #print(x_int) + + dx = abs(min([x_int - self.delta_stern(z), self.SB_profile(z) - x_int])) #Distance over which fillet occurs + + Xrad = [(x_int - dx*self.Rsb/2.0), ((x_int + dx*self.Rsb/2.0))] + + + # Parabolic Radius + + Yrad = [self.halfBeam_Stern([Xrad[0]], PROF)[0], self.halfBeam_SB(z, [Xrad[1]])[0]] + + dydx = [self.stern_dydx([Xrad[0]], PROF)[0], self.SB_dydx(z, [Xrad[1]])[0]] + + + + Arad = np.array([[Xrad[0]**4.0, Xrad[0]**3.0, Xrad[0]**2.0, Xrad[0], 1.0], + [Xrad[1]**4.0, Xrad[1]**3.0, Xrad[1]**2.0, Xrad[1], 1.0], + [4.0*Xrad[0]**3.0, 3.0*Xrad[0]**2.0, 2.0*Xrad[0], 1.0, 0.0], + [4.0*Xrad[1]**3.0, 3.0*Xrad[1]**2.0, 2.0*Xrad[1], 1.0, 0.0], + [12.0*Xrad[1]**2.0, 6.0*Xrad[1]**2.0, 2.0, 0.0, 0.0]]) + + + + brad = np.array([Yrad[0], Yrad[1], dydx[0], dydx[1], 0.0]) + + + PROFrad = np.linalg.solve(Arad,brad) + + + xstern = [i for i in x[:-1] if (i >= x1 and i <= Xrad[0])] + + xsbrad = [i for i in x[:-1] if (i > Xrad[0] and i < Xrad[1])] + + + xsb = [i for i in x[:-1] if (i >= Xrad[1] and i < x2)] + + ystern = self.halfBeam_Stern(xstern, PROF) + + + ysbrad = np.zeros((len(xsbrad),)) + + for i in range(0,len(xsbrad)): + + ysbrad[i] = PROFrad[0]*xsbrad[i]**4.0 + PROFrad[1]*xsbrad[i]**3.0 + PROFrad[2]*xsbrad[i]**2.0 + PROFrad[3]*xsbrad[i] + PROFrad[4] + + ysb = self.halfBeam_SB(z, xsb) + + XY[0:-1,0] = x[0:-1] + + XY[0:-1,1] = np.concatenate((ystern,ysbrad,ysb)) + + # Make sure XY[-1] is zero to create a closed mesh + XY[-1] = [x2,0.0] + + return XY + + def plot_BulbProfiles(self): + + # Plot intersection points in blue + # Plot chine pt in green + # Plot Center of Rc and Rk in red + # half Beam(z) in black + + z1 = np.linspace(0.0, self.WL, num = 200) + z2 = np.linspace(0.0, self.WL*self.HSBOA, num = 200) + x = np.zeros((200,2)) + for i in range(0,len(z1)): + x[i,0] = self.BB_profile(z1[i]) + x[i,1] = self.SB_profile(z2[i]) - self.Ls-self.Lm + + + fig2,ax2 = plt.subplots() + ax2.axis('equal') + #plt.axis([0,10,0,10]) + + ax2.plot(x[:,0], z1, '-', color = 'blue') + ax2.plot(x[:,1], z2,'-' ,color = 'red') + + + + + def BulbformConstraints(self): + + C = np.zeros((13,)) + + ''' + Bulb form constraints: + 0) BBlower z radius is less than Rk + 1) BB x radius is less than Rk + 2) BB half beam is less than the half beam of the midbody at the height of max BB beam + 3) BB is FWD of Delta Bow at 0 + 4) BB is FWD of Delta Bow at Vertex of Delta Bow + 5) BB is FWD of Delta Bow at WL + 6) SB lower z radius is less than Rk + 7) SB x radius is less than Rk + 8) SB half beam is less than the half beam of the midbody at the height of max BB beam + 9) SB Height overall is less than the height of the bottom of the transom + 10) delta_stern(z = 0) is fwd of the starting position of the max width of the bulb (SB_Prof[5]) + 11) delta_stern(z = WL*HSBOA) is fwd of the starting position of the max width of the bulb (SB_Prof[5]) + 12) if the vertex of delta_stern is less than WL*HSBOA, then it is also fwd of the starting position of the max width of the bulb (SB_Prof[5]) + ''' + + if self.bit_BB: + if self.Beta == 0.0: + C[0:3] = np.array([-1.0, + -1.0, + self.BB_Prof[3] - self.halfBeam_MidBody(self.BB_Prof[2])]) + + elif self.Rk > 0.0: + C[0:3] = np.array([self.BB_Prof[2] - self.Rk, + self.BB_Prof[3] -self.Rk, + self.BB_Prof[3] - self.halfBeam_MidBody(self.BB_Prof[2])]) + + else: + C[0:3] = np.array([1.0,1.0,1.0]) + + if self.DELTA_BOW[0] == 0.0: + Zv = -1.0 + else: + Zv = -self.DELTA_BOW[1]/ (2.0*self.DELTA_BOW[0]) + + if Zv >=0.0 and Zv <= self.WL: + vert_delta_bow = (self.delta_bow(Zv) - self.BB_Prof[5]) + else: + vert_delta_bow = -1.0 + + C[3:6] = np.array([self.BB_Prof[5] - self.delta_bow(0.0), + self.BB_Prof[5] - self.delta_bow(self.WL), + vert_delta_bow]) + + else: + C[0:6] = np.array([-1.0,-1.0,-1.0, -1.0, -1.0, -1.0]) + + + + if self.bit_SB: + if self.Beta == 0.0: + C[6:10] = np.array([-1.0, + -1.0, + self.SB_Prof[3] - self.halfBeam_MidBody(self.SB_Prof[2]), + self.WL*self.HSBOA - self.SK[1]]) + + elif self.Rk > 0.0: + C[6:10] = np.array([self.SB_Prof[2] - self.Rk, + self.SB_Prof[3] -self.Rk, + self.SB_Prof[3] - self.halfBeam_MidBody(self.SB_Prof[2]), + self.WL*self.HSBOA - self.SK[1]]) + + else: + C[6:10] = np.array([1.0,1.0,1.0,1.0]) + + if self.DELTA_STERN[0] == 0.0: + Zv = -1.0 + else: + Zv = -self.DELTA_STERN[1]/ (2.0*self.DELTA_STERN[0]) + + if Zv >=0.0 and Zv <= self.WL*self.HSBOA: + vert_delta_stern = (self.delta_stern(Zv) - self.SB_Prof[5]) + else: + vert_delta_stern = -1.0 + + C[10:13] = np.array([self.delta_stern(0.0) - self.SB_Prof[5], + self.delta_stern(self.WL*self.HSBOA) - self.SB_Prof[5], + vert_delta_stern]) + + + #If no Stern Bulb, then no constraint violations + else: + C[6:13] = np.array([-1.0,-1.0,-1.0, -1.0, -1.0, -1.0, -1.0]) + + + + + + + + return C + + + ''' + ===================================================================== + Section 6: Mesh Generatation + ====================================================================== + + This section contains functions needed to generate a complete mesh of the hull + as an STL + + + ''' + def gen_MeshGridPointCloud(self, NUM_WL = 51, PointsPerLOA = 501, Z = [], X = [], bit_GridOrList = 1): + # This generates each waterline with even x and z spacing in a grid + #Z and X assignments supercede NUM_WL and PointsPerLOA Assignments + + if len(Z) == 0: + Z = np.linspace(0.0001*self.Dd,self.Dd, NUM_WL) + + if len(X) == 0: + X = np.linspace(-self.LOA*0.5,1.5*self.LOA, 2*PointsPerLOA - 1) + + Points = [] + + for i in Z: + pts = self.gen_MeshGridWL(X,i) + + Points.append(pts) + + #returns the points if user wants a waterline structured array of points + if bit_GridOrList: + return Points + + #returns a list of points in shape (N,3) if a list is preferred + else: + Cloud = [] + for i in Points: + for j in i: + Cloud.append(j) + + return Cloud + + def gen_MeshGridWL(self, X, z): + + WL = [] + + if z == 0.0: + if self.bit_BB: + bow_start = self.BB_Prof[5] + else: + bow_start = self.Kappa_BOW*self.Lb + + if self.bit_SB: + stern_end = self.SB_Prof[5] + else: + stern_end = self.Lm + self.Lb + self.Kappa_STERN*self.Ls + + + WL.append([bow_start,0.0,0.0]) + + x = [i for i in X if (i > bow_start and i < stern_end)] + + for i in x: + WL.append([i,0.0,0.0]) + + WL.append([stern_end,0.0,0.0]) + + else: + # Now generate the remaining watelines + + # Figure out spacing and profiles + if self.bit_BB: + BOW = self.gen_waterline_bow_BB(z, X=X, bit_spaceOrGrid=0) + else: + BOW = self.gen_waterline_bow(z, X=X, bit_spaceOrGrid=0) + + for i in range(0,len(BOW)): + WL.append([BOW[i,0], BOW[i,1], z]) + + X_mid = [i for i in X if (i >= self.delta_bow(z) and i < self.delta_stern(z))] + Y_mid = self.halfBeam_MidBody(z) + + for i in range(0,len(X_mid)): + WL.append([X_mid[i], Y_mid, z]) + + if self.bit_SB: + STERN = self.gen_waterline_stern_SB(z, NUM_POINTS = 0, X=X, bit_spaceOrGrid=0) + else: + STERN = self.gen_waterline_stern(z, NUM_POINTS = 0, X=X, bit_spaceOrGrid=0) + + for i in range(0,len(STERN)): + WL.append([STERN[i,0], STERN[i,1], z]) + + return np.array(WL) + + + def gen_pointCloud(self, NUM_WL = 50, PointsPerWL = 300, bit_GridOrList = 0, Z = []): + # This function generates a point cloud [[xyz]0, ..., [xyz]N] of the hull + + #NUM_WL = number of waterlines to generate. The total waterlines must be greater than 8 to include these heights: z = [0, 0.001*self.Dd, BKz, SKz, Dd, WL, Hbb, and Hsb] + #PointsPerWL will be divided so that each major section of the ship (stern, midbody, and bow) each gets a proportion of the points to their total length + # bit_GridOrList determines of the PC will be returned as a grid (shape = (NUM_WL,PointsPerWL,3)) or a list (shape = (NUM_WL*PointsPerWL,3)) + # Z is an optional array of z heights to use for the point cloud or to generate the standard. Note that if Z is not empty, that the true number of waterlines will be equal to len(Z) + + if len(Z) == 0: + z = self.gen_WLHeights(NUM_WL) + else: + z = Z + NUM_WL = len(z) + + if bit_GridOrList: + PC = np.zeros((NUM_WL, PointsPerWL,3)) + + + for i in range(0,len(z)): + + PC[i] = self.gen_WLPoints(z[i], PointsPerWL) + + else: + PC = np.zeros((NUM_WL*PointsPerWL,3)) + for i in range(0,len(z)): + WL = self.gen_WLPoints(z[i], PointsPerWL) + PC[PointsPerWL*i:PointsPerWL*(i+1)] = np.array(WL) + + return PC + + def gen_WLHeights(self, NUM_WL): + + #NUM_WL = number of waterlines to generate. The total waterlines must be greater than 8 to include these heights: z = [0, 0.001*self.Dd, BKz, SKz, Dd, WL, Hbb, and Hsb] + + z = np.zeros((NUM_WL,)) + + z[0:7] = np.array([self.BK[1], self.SK[1], self.WL, self.Hbb*self.WL, self.HSBOA*self.WL, self.Hsb*self.WL*self.HSBOA, 0.001*self.Dd]) + z[7:] = np.linspace(0.0, self.Dd, NUM_WL - 7) + + z = np.sort(z) + + return z + + def gen_WLPoints(self, z, PointsPerWL = 300): + #This function generates the starboard waterline half breadths + #set up baseline first between deltabow and delta stern + # z is a scalar, PointsPerWL is also a scalar + + WL = [] + + if z == 0.0: + if self.bit_BB: + bow_start = self.BB_Prof[5] + else: + bow_start = self.Kappa_BOW*self.Lb + + if self.bit_SB: + stern_end = self.SB_Prof[5] + else: + stern_end = self.Lm + self.Lb + self.Kappa_STERN*self.Ls + + + + x = np.linspace(bow_start, stern_end, PointsPerWL) + for i in range(0,PointsPerWL): + WL.append([x[i], 0.0, 0.0]) + + else: + # Now generate the remaining watelines + if self.bit_BB and z <= self.WL: + bow_start = min([self.BB_profile(z), self.bow_profile(z)]) + else: + bow_start = self.bow_profile(z) + + if self.bit_SB and z <= self.HSBOA*self.WL: + + stern_end = max([self.SB_profile(z),self.stern_profile(z)]) + else: + stern_end = self.stern_profile(z) + + + WL_LOA = stern_end - bow_start + + #print([z,WL_LOA,stern_end, bow_start]) + + pts_bow = abs(int((self.delta_bow(z) - bow_start)/WL_LOA * PointsPerWL)) + 1 + pts_stern = abs(int((stern_end - self.delta_stern(z))/WL_LOA * PointsPerWL)) + 1 + + if (pts_bow + pts_stern) > PointsPerWL: + over = pts_bow + pts_stern - PointsPerWL + pts_mid = 0 + pts_stern = pts_stern - over + else: + pts_mid = PointsPerWL - pts_bow - pts_stern + + #print(self.delta_bow(z) - bow_start) + #print([z, pts_bow,pts_mid,pts_stern]) + + if self.bit_BB: + BOW = self.gen_waterline_bow_BB(z, NUM_POINTS = pts_bow) + else: + BOW = self.gen_waterline_bow(z, NUM_POINTS = pts_bow) + + for i in range(0,pts_bow): + WL.append([BOW[i,0], BOW[i,1], z]) + + X_mid = np.linspace(self.delta_bow(z), self.delta_stern(z), pts_mid+2) + Y_mid = self.halfBeam_MidBody(z) + + for i in range(0,pts_mid): + WL.append([X_mid[1+i], Y_mid, z]) + + if self.bit_SB: + STERN = self.gen_waterline_stern_SB(z, NUM_POINTS = pts_stern) + else: + STERN = self.gen_waterline_stern(z, NUM_POINTS = pts_stern) + + for i in range(0,pts_stern): + WL.append([STERN[i,0], STERN[i,1], z]) + + + + return np.array(WL) + + + def gen_stl(self, NUM_WL = 50, PointsPerWL = 300, bit_AddTransom = 1, bit_AddDeckLid = 0, bit_RefineBowAndStern = 0, namepath = 'Hull_Mesh') -> mesh.Mesh: + # This function generates a surface of the mesh by iterating through the points on the waterlines + + #compute number of triangles in the mesh + #hullTriangles = 2 * (2*PointsPerWL - 2) * (NUM_WL - 1) + #numTriangles = hullTriangles + transomTriangles = 0 + + #Generate WL + z = np.zeros((NUM_WL,)) + + z[0] = 0.0001*self.Dd + z[1] = 0.001*self.Dd + z[2:] = np.linspace(0.0, self.Dd, NUM_WL - 2) + + z = np.sort(z) + + x = np.linspace(-self.LOA*0.5,1.5*self.LOA, 2*PointsPerWL - 1) + + if bit_RefineBowAndStern: + # Add more points to X in the bow and stern + + x_sub1 = x[0:int(0.75*PointsPerWL)] + 0.5*(x[1] - x[0]) + x_sub2 = x[-int(0.75*PointsPerWL):] + 0.5*(x[1] - x[0]) + x = np.concatenate((x_sub1, x_sub2, x)) + x = np.sort(x) + + + + + #Generate MeshGrid PC + pts = self.gen_MeshGridPointCloud(NUM_WL = NUM_WL, PointsPerLOA = PointsPerWL, Z = z, X = x, bit_GridOrList = 1) + + #start to assemble the triangles into vectors of indices from pts + TriVec = [] + + for i in range(0,NUM_WL-1): + + #Find idx where the mesh grids begin to align between two rows returns a zero or 1: + + bow = np.argmax([pts[i][0,0],pts[i+1][0,0]]) + + stern = np.argmin([pts[i][-1,0],pts[i+1][-1,0]]) + + + + # Find index where mesh grid lines up and ends between each WL + + if bow: + idx_WLB1 = 1 + idx_WLB0 = np.where(pts[i][:,0] == pts[i+1][idx_WLB1,0])[0][0] + else: + idx_WLB0 = 1 + idx_WLB1 = np.where(pts[i+1][:,0] == pts[i][idx_WLB0,0])[0][0] + + if stern: + idx_WLS1 = len(pts[i+1]) - 2 + idx_WLS0 = np.where(pts[i][:,0] == pts[i+1][idx_WLS1,0])[0][0] + else: + idx_WLS0 = len(pts[i]) - 2 + idx_WLS1 = np.where(pts[i+1][:,0] == pts[i][idx_WLS0,0])[0][0] + + #check that these two are the same size: + + #Build the bow triangles Includes Port assignments + + if bow: + TriVec.append([pts[i+1][idx_WLB1], pts[i][0], pts[i+1][0]]) + + for j in range(0,idx_WLB0): + TriVec.append([pts[i+1][idx_WLB1], pts[i][j+1], pts[i][j]]) + + + + else: + + for j in range(0,idx_WLB1): + TriVec.append([pts[i][0],pts[i+1][j], pts[i+1][j+1]]) + + TriVec.append([pts[i][0],pts[i+1][idx_WLB1], pts[i][idx_WLB0]]) + + #Build main part of hull triangles. Port Assignments + for j in range(0, idx_WLS1-idx_WLB1): + + TriVec.append([pts[i][idx_WLB0+j], pts[i+1][idx_WLB1+j], pts[i+1][idx_WLB1+j+1]]) + TriVec.append([pts[i][idx_WLB0+j], pts[i+1][idx_WLB1+j+1], pts[i][idx_WLB0+j+1]]) + + #Build the stern: + if stern: + + for j in range(idx_WLS0,len(pts[i])-1): + TriVec.append([pts[i+1][idx_WLS1], pts[i][j+1],pts[i][j]]) + + TriVec.append([pts[i+1][idx_WLS1], pts[i+1][-1], pts[i][-1]]) + + else: + + TriVec.append([pts[i][idx_WLS0], pts[i+1][idx_WLS1], pts[i][-1]]) + + for j in range(idx_WLS1, len(pts[i+1])-1): + TriVec.append([pts[i][-1], pts[i+1][j], pts[i+1][j+1]]) + + + TriVec = np.array(TriVec) + + hullTriangles = 2*len(TriVec) + numTriangles = hullTriangles + + + + #add triangles if there is a transom + if bit_AddTransom: + wl_above = len([i for i in z if i > self.SK[1]]) + + z_idx = NUM_WL - wl_above - 1 + + transomTriangles = 2*wl_above - 1 + + numTriangles += transomTriangles + + #Add triangles if there is a deck lid (meaning the ship becomes a closed body) + if bit_AddDeckLid: + numTriangles += 2*len(pts[-1]) - 3 + + + HULL = mesh.Mesh(np.zeros(numTriangles, dtype=mesh.Mesh.dtype)) + + HULL.vectors[0:len(TriVec)] = np.copy(TriVec) + + TriVec_stbd = np.copy(TriVec[:,::-1]) + TriVec_stbd[:,:,1] *= -1 + HULL.vectors[len(TriVec):hullTriangles] = np.copy(TriVec_stbd) + + # NowBuild the transom: + if bit_AddTransom: + + + pts_trans = np.zeros((wl_above+1,3)) + + for i in range(0,len(pts_trans)): + pts_trans[i] = pts[z_idx+i][-1,:] + + + + pts_tranp = np.array(pts_trans) + + pts_tranp[:,1] *= -1.0 + + + + + HULL.vectors[hullTriangles] = np.array([pts_trans[0], pts_trans[1], pts_tranp[1]]) + + for i in range(1,wl_above): + HULL.vectors[hullTriangles + 2*i-1] = np.array([pts_trans[i], pts_trans[i+1], pts_tranp[i]]) + HULL.vectors[hullTriangles + 2*i] = np.array([pts_tranp[i], pts_trans[i+1], pts_tranp[i+1]]) + + + # Add the deck lid + if bit_AddDeckLid: + + #pts_Lids are starboard points on the deck + #pts_Lidp are port points on the deck + + pts_Lids = pts[NUM_WL-1] + + pts_Lidp = np.array(pts_Lids) + pts_Lidp[:,1] *= -1.0 + + startTriangles = hullTriangles + transomTriangles + + # Points are orered so the right hand rule points the lid in positive z + HULL.vectors[startTriangles] = np.array([pts_Lids[0], pts_Lidp[1], pts_Lids[1]]) + + for i in range(1,len(pts_Lids)-1): + HULL.vectors[startTriangles + 2*i - 1] = np.array([pts_Lids[i], pts_Lidp[i], pts_Lids[i+1]]) + HULL.vectors[startTriangles + 2*i] = np.array([pts_Lids[i+1], pts_Lidp[i], pts_Lidp[i+1]]) + + if not (namepath is None): + HULL.save(namepath + '.stl') + return HULL + + + + def gen_PC_for_Cw(self, draft, NUM_WL = 51, PointsPerWL = 301): + ''' + This code generates the Point Grid and the Inputs for the Cw prediction. + 0) Z and X -> Translated X and Z vectors for the points used in the solver. + X[0] = 0, X[-1] = LOA of submerged Hull. Z[0] = -draft Z[-1] = 0 + 1) WL = length of Waterline + 2) Y -> point grid of Y offsets + + ''' + Z = np.linspace(0.00000001*self.Dd, draft, NUM_WL) + + x_bow = np.zeros((len(Z),)) + x_stern = np.zeros((len(Z),)) + + for i in range(0,len(Z)): + # Now generate the remaining watelines + if self.bit_BB and Z[i] <= self.WL: + x_bow[i] = self.BB_profile(Z[i]) + else: + x_bow[i] = self.bow_profile(Z[i]) + + if self.bit_SB and Z[i] <= self.HSBOA*self.WL: + + x_stern[i] = self.SB_profile(Z[i]) + else: + x_stern[i] = self.stern_profile(Z[i]) + + + WL = x_stern[-1] - x_bow[-1] + + X = np.linspace(np.amin(x_bow), np.amax(x_stern), PointsPerWL) + + Y = np.zeros((PointsPerWL,NUM_WL)) + + points = self.gen_MeshGridPointCloud(Z = Z, X = X, bit_GridOrList = 1) + + + + for i in range(0,len(Z)): + idx = np.where(X == points[i][1][0])[0][0] #points[i,1,0] = first X in points where y != 0 + + + + for j in range(1,len(points[i])-1): + Y[idx+j-1,i] = points[i][j][1] + + + X = X - X[0] # Normalize so that X[0] = 0 + Z = Z - Z[-1] # Normalize so that Z[-1] = 0 + + return X,Z,Y,WL + + + + + + + + def input_Constraints(self): + + return np.concatenate((self.GenralHullformConstraints(),self.CrossSectionConstraints(), self.BowformConstraints(), self.SternformConstraints(), self.BulbformConstraints())) + + + + ''' + ========================================================================= + Section 7: Geometric and Volumetric Analysis Fucntions + ========================================================================== + + This section contains functions to perform Naval Architecture related analyis on the geometry of the hull + + 0) Displacement(z) -> calculates the submerged volume at a given height, z + 1) CentersOfBuoyancy(z) -> Returns the centers of Buoyancy for a given waterline height, z + 1) WaterplaneArea(z) -> caluclates the area of the waterplane at height, z + 2) WaterplaneMoments(z) -> calculates the Center of Flotation, Ixx, and Iyy of the Waterplane at height(z) + 3) MTC(z) -> Not Implemented yet + 4) Righting Moment(z) -> Not Implemented Yet + 5) Block Coefficient -> Not Implemented yet + 6) Draft, Heel and Trim -> Not Implemented Yet + 7) LOA_wBulb -> Returns the maximum length of the hull including added lengths from bulbs + 8) Max_Beam_midship -> Returns the maximum beam of the midship section (calculated from midship section functions) + 9) Max_Beam_PC -> Returns the maximum beam of the hull (Estimated from point cloud for volume calculations) + + ''' + + def Calc_VolumeProperties(self, NUM_WL = 101, PointsPerWL = 1000): + #This function generates a point cloud to be used for volumetric measurements an calls all of the evaluation measurements to be calculated. + Z = np.linspace(0.00001,self.Dd, num=NUM_WL) + + self.PCMeasurement = self.gen_pointCloud(NUM_WL = NUM_WL, PointsPerWL = PointsPerWL, bit_GridOrList = 1, Z = Z) + + self.Calc_WaterPlaneArea() + + self.Calc_Volumes() + + self.Calc_LCFs() + + self.Calc_CB(Z) + + self.Calc_2ndMoments() + + self.Calc_WettedSurface(Z) + + self.Calc_WaterlineLength() + + return Z + + + def Calc_Volumes(self): + # this Function calculates the 3D Volume of a hull below a height z by integrating the waterplane areas of each waterline below z as well. + Vol = np.zeros((len(self.PCMeasurement),)) + + Vol[0] = 0.5*self.Areas_WP[0]*self.PCMeasurement[0,0,2] + + for i in range(1,len(Vol)): + Vol[i] = Vol[i-1] + 0.5*(self.Areas_WP[i] + self.Areas_WP[i-1])*(self.PCMeasurement[i,0,2] - self.PCMeasurement[i-1,0,2]) + + self.Volumes = Vol + + def Calc_WaterPlaneArea(self): + # This function calcuates the waterplane area for a given z height + Areas = np.zeros((len(self.PCMeasurement),)) + + for i in range(0,len(Areas)): + + WL = self.PCMeasurement[i] + Areas[i] = 2.0*np.trapz(WL[:,1], x=WL[:,0]) + + self.Areas_WP = Areas + + def Calc_LCFs(self): + #This function calculates the Longitudinal Center of Flotation for the Waterplane sections + LCF = np.zeros((len(self.PCMeasurement),)) + + for i in range(0,len(LCF)): + + Moment = 0 + + for j in range(1,len(self.PCMeasurement[i])): + #sum up trapezoid moments of area + + + dx = self.PCMeasurement[i,j,0] - self.PCMeasurement[i,j-1,0] + a = 2.0 * self.PCMeasurement[i,j,1] + b = 2.0 * self.PCMeasurement[i,j-1,1] + + + + cx = self.PCMeasurement[i,j-1,0] + dx/3.0 * (2.0*a + b) / (a + b) + + Moment = Moment + cx * 0.5*(a+b)*dx + + LCF[i] = Moment/self.Areas_WP[i] + + self.LCFs = LCF + + def Calc_CB(self,Z): + #This function calculates the longitudinal and vertical centers of buoyancy for each of the waterlines + #CB[:,0] provides the LCB and CB[:,1] is the VCB for eah volume + + + CB = np.zeros((len(self.PCMeasurement),2)) + + #Calculate Moments for X and Z directions for the Center of Buoynacy + + MomentX = np.multiply(self.LCFs, self.Areas_WP) + + MomentZ = np.multiply(Z,self.Areas_WP) + + + for i in range(0,len(CB)): + + CB[i,0] = np.trapz(MomentX[0:i+1],x = Z[0:i+1])/self.Volumes[i] + CB[i,1] = np.trapz(MomentZ[0:i+1],x = Z[0:i+1])/self.Volumes[i] + + self.VolumeCentroids = CB + + + def Calc_2ndMoments(self): + #this function calculates the second moment of area Ixx and Iyy for each WaterPlane in the form I[i] = [Ixxi,Iyyi] + + I = np.zeros((len(self.PCMeasurement),2)) + + for i in range(0,len(I)): + + Ixx = 0.0 + Iyy = 0.0 + + for j in range(1,len(self.PCMeasurement[i])): + #sum up trapezoid moments of area + dx = self.PCMeasurement[i,j,0] - self.PCMeasurement[i,j-1,0] + a = 2.0 * self.PCMeasurement[i,j,1] + b = 2.0 * self.PCMeasurement[i,j-1,1] + + cx = self.PCMeasurement[i,j-1,0] + dx/3.0 * (2.0*a+b)/(a+b) + + Ixx = Ixx + dx/48.0 * (a + b) * (a**2.0 + b**2.0) + + Iyy = Iyy + dx**3.0 * (a**2.0 + 4.0*a*b + b**2.0)/(36*(a+b)) + 0.5*(a+b)*dx*(self.LCFs[i] - cx)**2.0 + + I[i] = [Ixx,Iyy] + + self.I_WP = I + + def Calc_WettedSurface(self,Z): + # This function calcultes and summates the wetted surface between each draft line (Z[]and at the bottom of the hull, by estimating length along the surface + + ArcL = np.zeros((len(Z),)) + WSA = np.zeros((len(Z),)) + + + for i in range(0,len(ArcL)): + + for j in range(1,len(self.PCMeasurement[0])): + + #dL = distance of length along outside of ship along waterline at z[i] + dL = np.sqrt((self.PCMeasurement[i,j,0] - self.PCMeasurement[i,j-1,0])**2.0 + (self.PCMeasurement[i,j,1] - self.PCMeasurement[i,j-1,1])**2.0) + + ArcL[i] = ArcL[i] + dL + + ArcL[i] + ArcL[i] + self.PCMeasurement[i,-1,1] #Add transom width to Arc Length + + #Wetted Surface Area is integral of Arc Length from 0 to Z[idx] + + WSA[0] = 2*self.Areas_WP[0] #wetted surface area at bottom of hull is approximately area of waterplane at z ~=0 + + for i in range(1,len(Z)): + + #Wetted Surface area is cumulative sum of WSA at each height (2*0.5 for trapezoid rule and for 2 sides of hull = 1) + WSA[i] = WSA[i-1] + (ArcL[i]+ArcL[i-1])*(Z[i] - Z[i-1]) + + self.Area_WS = WSA + + + def Calc_WaterlineLength(self): + #This function returns the length of the waterline for each Z + + WLL = np.zeros((len(self.PCMeasurement),)) + + for i in range(0,len(WLL)): + #Length of Stern Position - Bow Position in X + WLL[i] = self.PCMeasurement[i,-1,0] - self.PCMeasurement[i,0,0] + + self.WL_Lengths = WLL + + def Calc_LOA_wBulb(self): + #This function returns the length of the hull including the bulb lengths + + if self.bit_BB: + bow_start = min([0.0,self.BB_Prof[5]-self.BB_Prof[4]]) + else: + bow_start = 0.0 + + if self.bit_SB: + stern_end = max([self.LOA, self.SB_Prof[5]+self.SB_Prof[4]]) + else: + stern_end = self.LOA + + self.LOA_wBulb = stern_end - bow_start + + return self.LOA_wBulb + + def Calc_Max_Beam_midship(self): + #This function returns the maximum beam of the midship section (calculated from midship section functions) + + #fist check Bd vs Bc + if self.Bd >= self.Bc: + self.Max_Beam_midship = self.Bd*2.0 + else: + self.Max_Beam_midship = (self.Rc_Center[0] + self.Rc)*2.0 #Max beam is the y coordinate of the center of the chine plus the radius of the chine + + return self.Max_Beam_midship + + def Calc_Max_Beam_PC(self): + #This function returns the maximum beam of the hull (Estimated from point cloud for volume calculations) + + self.Max_Beam_PC = 2.0*np.amax(self.PCMeasurement[:,:,1]) + + return self.Max_Beam_PC + + + def interp(A,Z,z): + # This function interpolates data to approximate A(z) given values of A(Z) + + idx = np.where(Z < z)[0][-1] + + frac = (z - Z[idx])/(Z[idx+1] - Z[idx]) + + return A[idx] + frac*(A[idx+1] - A[idx]) + +