In [4]:
import numpy as np

%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt

mpl.rcParams['figure.figsize'] = (16, 9)

import sys
import imp
import os

import polymesh.mesh as mesh
import polymesh.hydrostatic as hydrostatic

In [15]:
# Import the test geoemtry, which is a full size version of the ship KCS
shipMesh = mesh.importObj('BrAa_40m_hull_side_offset.obj')

# The keel is located at z=0 in the geometry file. Translate it downwards to right depth
T = 1.65
shipMesh.translate(0, 0, -T)

In [16]:
rho = 1025

wetMesh = hydrostatic.extractWetSurface(shipMesh)
wetMesh.exportObj('wetMesh.obj')
wetMesh.calculateEdgeConnectivity()

Volume         = hydrostatic.calculateVolume(wetMesh)
Surface        = hydrostatic.calculateSurface(wetMesh)
volumeCentroid = hydrostatic.calculateVolumeCentroid(wetMesh)
Dimensions     = hydrostatic.calculateDimensions(wetMesh) 
Mass = Volume*rho

print('Volume:', np.round(Volume, decimals=3), 'm^3')
print('Surface:', np.round(Surface, decimals=3), 'm^2')
print('Volume centroid:', np.round(volumeCentroid[0], decimals=4), np.round(volumeCentroid[1], decimals=4), np.round(volumeCentroid[2], decimals=4))
print('Dimensions:', Dimensions)
print('Mass:', Mass/1e3, 'tonnes')


Volume: 102.642 m^3
Surface: 173.284 m^2
Volume centroid: 16.2849 3.6486 -0.682
Dimensions: [40.        2.644972  0.640746]
Mass: 105.20835872925721 tonnes

In [ ]:
shipMesh.translate(0, 0, T)

In [62]:
n_planes    = 51
norm_axis   = 0
sort_axis   = 1

mesh_       = wetMesh

plt.scatter(mesh_.verts[mesh_.face_verts][:,0], mesh_.verts[mesh_.face_verts][:,1])
ax_min      = np.amin(mesh_.verts[mesh_.face_verts][:,0])
ax_max      = np.amax(mesh_.verts[mesh_.face_verts][:,0])
ax_len      = ax_max-ax_min
d_long      = (ax_len)/(n_planes-1)
pos_long    = np.zeros(n_planes)
plane_points= []


for i in range(n_planes):
    pos_long[i] = ax_min + i*d_long
    plane       = hydrostatic.extractPlane(wetMesh, norm_axis, pos_long[i])
    points      = plane[0]
    points      = points[points[:,sort_axis].argsort()]
    plane_points.append(points)
    
    
print(plane_points[0])


[[0.       2.327585 1.65    ]
 [0.       2.338623 1.545334]
 [0.       2.349707 1.440971]
 [0.       2.360959 1.337098]
 [0.       2.372801 1.233434]
 [0.       2.38587  1.12948 ]
 [0.       2.401368 1.02433 ]
 [0.       2.422419 0.918771]
 [0.       2.459284 0.820772]
 [0.       2.521062 0.738329]
 [0.       2.605076 0.672343]
 [0.       2.704853 0.621661]
 [0.       2.810518 0.584013]
 [0.       2.913747 0.556838]
 [0.       3.015817 0.537549]
 [0.       3.119643 0.523832]
 [0.       3.225114 0.514498]
 [0.       3.331438 0.508503]
 [0.       3.438149 0.504247]
 [0.       3.543825 0.503613]
 [0.       3.6495   0.50298 ]
 [0.       3.64975  0.50298 ]
 [0.       3.65     0.50298 ]
 [0.       3.65025  0.50298 ]
 [0.       3.6505   0.50298 ]
 [0.       3.756175 0.503613]
 [0.       3.861851 0.504247]
 [0.       3.968562 0.508503]
 [0.       4.074886 0.514498]
 [0.       4.180357 0.523832]
 [0.       4.284184 0.537549]
 [0.       4.386253 0.556838]
 [0.       4.489482 0.584013]
 [0.       4.595148 0.621661]
 [0.       4.694924 0.672343]
 [0.       4.778938 0.738329]
 [0.       4.840716 0.820772]
 [0.       4.877581 0.918771]
 [0.       4.898632 1.02433 ]
 [0.       4.91413  1.12948 ]
 [0.       4.927199 1.233434]
 [0.       4.939041 1.337098]
 [0.       4.950293 1.440971]
 [0.       4.961378 1.545334]
 [0.       4.972415 1.65    ]]

In [67]:
plane_no = 49
plt.figure()
plt.plot(plane_points[plane_no][:,1],plane_points[plane_no][:,2])
plt.tight_layout()



In [69]:
vessel_name = 'BrAa_40m_hull'
f           = open(vessel_name+'.mgf', 'w')
f.write(vessel_name+'\n')
for i in range(3):
    f.write('\n')
f.write(str(ax_len)+'\n')

for i in range(n_planes-1):
    plane_index     = n_planes-2-i
    shpx_index      = plane_index+1
    
    f.write(str(shpx_index)+'\n')
    f.write(str(pos_long[plane_index]-.5*ax_len)+'\n')
    n_points_plane = len(plane_points[plane_index])
    f.write(str(n_points_plane)+'\n')
    for j in range(n_points_plane):
        point_index     = n_points_plane-1-j
        f.write('{:.8f}    {:.8f}\n'.format(plane_points[plane_index][point_index,1], plane_points[plane_index][point_index,2]))
    


f.close()

In [ ]: