From 9770cca561854c7f9d0d7e75b29863f9956fd88b Mon Sep 17 00:00:00 2001 From: Jose Luis Cercos Pita Date: Wed, 20 Jan 2016 15:13:00 +0100 Subject: [PATCH] Finished the GZ curve implementation --- src/Mod/Ship/TankInstance.py | 103 +++++++++++++++++++++-- src/Mod/Ship/shipGZ/PlotAux.py | 139 +++++++++++-------------------- src/Mod/Ship/shipGZ/TaskPanel.py | 12 +-- src/Mod/Ship/shipGZ/Tools.py | 37 +++++--- 4 files changed, 178 insertions(+), 113 deletions(-) diff --git a/src/Mod/Ship/TankInstance.py b/src/Mod/Ship/TankInstance.py index aa2dafcd0..038e6e918 100644 --- a/src/Mod/Ship/TankInstance.py +++ b/src/Mod/Ship/TankInstance.py @@ -26,7 +26,7 @@ from math import * from PySide import QtGui, QtCore import FreeCAD as App import FreeCADGui as Gui -from FreeCAD import Base, Vector, Placement, Rotation +from FreeCAD import Base, Vector, Matrix, Placement, Rotation import Part import Units from shipUtils import Paths, Math @@ -79,14 +79,23 @@ class Tank: """ pass - def getVolume(self, fp, level): + def getVolume(self, fp, level, return_shape=False): """Return the fluid volume inside the tank, provided the filling level. Keyword arguments: fp -- Part::FeaturePython object affected. level -- Percentage of filling level (interval [0, 1]). + return_shape -- False if the tool should return the fluid volume value, + True if the tool should return the volume shape. """ - max(min(level, 1.0), 0.0) + if level <= 0.0: + if return_shape: + return Part.Vertex() + return Units.Quantity(0.0, Units.Volume) + if level >= 1.0: + if return_shape: + return fp.Shape.copy() + return Units.Quantity(fp.Shape.Volume, Units.Volume) # Build up the cutting box bbox = fp.Shape.BoundBox @@ -97,7 +106,7 @@ class Tank: box = App.ActiveDocument.addObject("Part::Box","Box") length_format = USys.getLengthFormat() box.Placement = Placement(Vector(bbox.XMin - dx, - bbox.XMin - dy, + bbox.YMin - dy, bbox.ZMin - dz), Rotation(App.Vector(0,0,1),0)) box.Length = length_format.format(3.0 * dx) @@ -133,7 +142,91 @@ class Tank: random_bounds)) App.ActiveDocument.recompute() - return Units.Quantity(common.Shape.Volume, Units.Volume) + if return_shape: + ret_value = common.Shape.copy() + else: + ret_value = Units.Quantity(common.Shape.Volume, Units.Volume) + + App.ActiveDocument.removeObject(common.Name) + App.ActiveDocument.removeObject(tank.Name) + App.ActiveDocument.removeObject(box.Name) + App.ActiveDocument.recompute() + + return ret_value + + def getCoG(self, fp, vol, roll=0.0, trim=0.0): + """Return the fluid volume center of gravity, provided the volume of + fluid inside the tank. + + The returned center of gravity is refered to the untransformed ship. + + Keyword arguments: + fp -- Part::FeaturePython object affected. + vol -- Volume of fluid (in m^3). + roll -- Ship roll angle (in degrees). + trim -- Ship trim angle (in degrees). + + If the fluid volume is bigger than the total tank one, it will be + conveniently clamped. + """ + # Change the units of the volume, and clamp the value + vol = vol * Units.Metre.Value**3 + if vol <= 0.0: + return Vector() + if vol >= fp.Shape.Volume: + vol = 0.0 + for solid in fp.Shape.Solids: + vol += solid.Volume + sCoG = solid.CenterOfMass + cog.x = cog.x + sCoG.x * solid.Volume + cog.y = cog.y + sCoG.y * solid.Volume + cog.z = cog.z + sCoG.z * solid.Volume + cog.x = cog.x / vol + cog.y = cog.y / vol + cog.z = cog.z / vol + return cog + + # Get a first estimation of the level + level = vol / fp.Shape.Volume + + # Transform the tank shape + current_placement = fp.Placement + m = current_placement.toMatrix() + m.rotateX(radians(roll)) + m.rotateY(-radians(trim)) + fp.Placement = m + + # Iterate to find the fluid shape + for i in range(COMMON_BOOLEAN_ITERATIONS): + shape = self.getVolume(fp, level, return_shape=True) + error = (vol - shape.Volume) / fp.Shape.Volume + if abs(error) < 0.01: + break + level += error + + # Get the center of gravity + vol = 0.0 + cog = Vector() + if len(shape.Solids) > 0: + for solid in shape.Solids: + vol += solid.Volume + sCoG = solid.CenterOfMass + cog.x = cog.x + sCoG.x * solid.Volume + cog.y = cog.y + sCoG.y * solid.Volume + cog.z = cog.z + sCoG.z * solid.Volume + cog.x = cog.x / vol + cog.y = cog.y / vol + cog.z = cog.z / vol + + # Untransform the object to retrieve the original position + fp.Placement = current_placement + p = Part.Point(cog) + m = Matrix() + m.rotateY(radians(trim)) + m.rotateX(-radians(roll)) + p.rotate(Placement(m)) + + return Vector(p.X, p.Y, p.Z) class ViewProviderTank: diff --git a/src/Mod/Ship/shipGZ/PlotAux.py b/src/Mod/Ship/shipGZ/PlotAux.py index 4d6427a71..cbbf56879 100644 --- a/src/Mod/Ship/shipGZ/PlotAux.py +++ b/src/Mod/Ship/shipGZ/PlotAux.py @@ -22,36 +22,37 @@ #*************************************************************************** import os +import math from PySide import QtGui, QtCore import FreeCAD import FreeCADGui -from FreeCAD import Base import Spreadsheet +from shipUtils import Paths + class Plot(object): - def __init__(self, x, y, disp, xcb, ship): - """ Constructor. performs the plot and shows it. - @param x X coordinates. - @param y Transversal computed areas. - @param disp Ship displacement. - @param xcb Bouyancy center length. - @param ship Active ship instance. - """ - self.plot(x, y, disp, xcb, ship) - self.spreadSheet(x, y, ship) + def __init__(self, roll, gz, draft, trim): + """ Plot the GZ curve - def plot(self, x, y, disp, xcb, ship): - """ Perform the areas curve plot. - @param x X coordinates. - @param y Transversal areas. - @param disp Ship displacement. - @param xcb Bouyancy center length. - @param ship Active ship instance. - @return True if error happens. + Position arguments: + roll -- List of roll angles (in degrees). + gz -- List of GZ values (in meters). + draft -- List of equilibrium drafts (in meters). + trim -- List of equilibrium trim angles (in degrees). + """ + self.plot(roll, gz) + self.spreadSheet(roll, gz, draft, trim) + + def plot(self, roll, gz): + """ Plot the GZ curve. + + Position arguments: + roll -- List of roll angles (in degrees). + gz -- List of GZ values (in meters). """ try: import Plot - plt = Plot.figure('Areas curve') + plt = Plot.figure('GZ') except ImportError: msg = QtGui.QApplication.translate( "ship_console", @@ -60,90 +61,46 @@ class Plot(object): QtGui.QApplication.UnicodeUTF8) FreeCAD.Console.PrintWarning(msg + '\n') return True - # Plot areas curve - areas = Plot.plot(x, y, 'Transversal areas') - areas.line.set_linestyle('-') - areas.line.set_linewidth(2.0) - areas.line.set_color((0.0, 0.0, 0.0)) - # Get perpendiculars data - Lpp = ship.Length.getValueAs('m').Value - FPx = 0.5 * Lpp - APx = -0.5 * Lpp - maxArea = max(y) - # Plot perpendiculars - FP = Plot.plot([FPx, FPx], [0.0, maxArea]) - FP.line.set_linestyle('-') - FP.line.set_linewidth(1.0) - FP.line.set_color((0.0, 0.0, 0.0)) - AP = Plot.plot([APx, APx], [0.0, maxArea]) - AP.line.set_linestyle('-') - AP.line.set_linewidth(1.0) - AP.line.set_color((0.0, 0.0, 0.0)) - # Add annotations for prependiculars + + gz_plot = Plot.plot(roll, gz, 'GZ curve') + gz_plot.line.set_linestyle('-') + gz_plot.line.set_linewidth(1.0) + gz_plot.line.set_color((0.0, 0.0, 0.0)) + ax = Plot.axes() - ax.annotate('AP', xy=(APx + 0.01 * Lpp, 0.01 * maxArea), size=15) - ax.annotate('AP', xy=(APx + 0.01 * Lpp, 0.95 * maxArea), size=15) - ax.annotate('FP', xy=(FPx + 0.01 * Lpp, 0.01 * maxArea), size=15) - ax.annotate('FP', xy=(FPx + 0.01 * Lpp, 0.95 * maxArea), size=15) - # Add some additional data - addInfo = ("$XCB = {0} \\; \\mathrm{{m}}$\n" - "$Area_{{max}} = {1} \\; \\mathrm{{m}}^2$\n" - "$\\bigtriangleup = {2} \\; \\mathrm{{tons}}$".format( - xcb, - maxArea, - disp)) - ax.text(0.0, - 0.01 * maxArea, - addInfo, - verticalalignment='bottom', - horizontalalignment='center', - fontsize=20) - # Write axes titles - Plot.xlabel(r'$x \; \mathrm{m}$') - Plot.ylabel(r'$Area \; \mathrm{m}^2$') + Plot.xlabel(r'$\phi \; [\mathrm{deg}]$') + Plot.ylabel(r'$GZ \; [\mathrm{m}]$') ax.xaxis.label.set_fontsize(20) ax.yaxis.label.set_fontsize(20) - # Show grid + Plot.grid(True) - # End plt.update() return False - def spreadSheet(self, x, y, ship): - """ Write the output data file. - @param x X coordinates. - @param y Transversal areas. - @param ship Active ship instance. + def spreadSheet(self, roll, gz, draft, trim): + """ Create a Spreadsheet with the results + + Position arguments: + roll -- List of roll angles (in degrees). + gz -- List of GZ values (in meters). + draft -- List of equilibrium drafts (in meters). + trim -- List of equilibrium trim angles (in degrees). """ s = FreeCAD.activeDocument().addObject('Spreadsheet::Sheet', - 'Areas curve') + 'GZ') # Print the header - s.set("A1", "x [m]") - s.set("B1", "area [m^2]") - s.set("C1", "FP x") - s.set("D1", "FP y") - s.set("E1", "AP x") - s.set("F1", "AP y") + s.set("A1", "roll [deg]") + s.set("B1", "GZ [m]") + s.set("C1", "draft [m]") + s.set("D1", "trim [deg]") - # Print the perpendiculars data - Lpp = ship.Length.getValueAs('m').Value - FPx = 0.5 * Lpp - APx = -0.5 * Lpp - maxArea = max(y) - s.set("C2", str(FPx)) - s.set("D2", str(0.0)) - s.set("C3", str(FPx)) - s.set("D3", str(maxArea)) - s.set("E2", str(APx)) - s.set("F2", str(0.0)) - s.set("E3", str(APx)) - s.set("F3", str(maxArea)) - # Print the data - for i in range(len(x)): - s.set("A{}".format(i + 2), str(x[i])) - s.set("B{}".format(i + 2), str(y[i])) + for i in range(len(roll)): + s.set("A{}".format(i + 2), str(roll[i])) + s.set("B{}".format(i + 2), str(gz[i])) + s.set("C{}".format(i + 2), str(draft[i])) + s.set("D{}".format(i + 2), str(trim[i])) # Recompute FreeCAD.activeDocument().recompute() \ No newline at end of file diff --git a/src/Mod/Ship/shipGZ/TaskPanel.py b/src/Mod/Ship/shipGZ/TaskPanel.py index 496b5b516..64173cd97 100644 --- a/src/Mod/Ship/shipGZ/TaskPanel.py +++ b/src/Mod/Ship/shipGZ/TaskPanel.py @@ -55,11 +55,13 @@ class TaskPanel: for i in range(n_points): rolls.append(roll * i / float(n_points - 1)) - gz = Tools.solve(self.ship, - self.weights, - self.tanks, - rolls, - form.var_trim.isChecked()) + gzs, drafts, trims = Tools.solve(self.ship, + self.weights, + self.tanks, + rolls, + form.var_trim.isChecked()) + + PlotAux.Plot(rolls, gzs, drafts, trims) return True diff --git a/src/Mod/Ship/shipGZ/Tools.py b/src/Mod/Ship/shipGZ/Tools.py index c48514551..cf5d3b230 100644 --- a/src/Mod/Ship/shipGZ/Tools.py +++ b/src/Mod/Ship/shipGZ/Tools.py @@ -47,7 +47,8 @@ def solve(ship, weights, tanks, rolls, var_trim=True): @param rolls List of considered roll angles. @param var_trim True if the trim angle should be recomputed at each roll angle, False otherwise. - @return GZ values for each roll angle + @return GZ values, drafts and trim angles, for each roll angle (in 3 + separated lists) """ # Get the unloaded weight (ignoring the tanks for the moment). W = 0.0 @@ -63,33 +64,43 @@ def solve(ship, weights, tanks, rolls, var_trim=True): # Get the tanks weight TW = 0.0 + VOLS = [] for t in tanks: # t[0] = tank object # t[1] = load density # t[2] = filling level vol = t[0].Proxy.getVolume(t[0], t[2]).getValueAs('m^3').Value + VOLS.append(vol) TW += vol * t[1] TW = TW * G gzs = [] + drafts = [] + trims = [] for i,roll in enumerate(rolls): App.Console.PrintMessage("{0} / {1}\n".format(i + 1, len(rolls))) - gz = solve_point(W, COG, TW, ship, tanks, roll, var_trim) + gz, draft, trim = solve_point(W, COG, TW, VOLS, + ship, tanks, roll, var_trim) if gz is None: return [] gzs.append(gz) + drafts.append(draft) + trims.append(trim) - return gzs + return gzs, drafts, trims -def solve_point(W, COG, TW, ship, tanks, roll, var_trim=True): +def solve_point(W, COG, TW, VOLS, ship, tanks, roll, var_trim=True): """ Compute the ship GZ value. @param W Empty ship weight. @param COG Empty ship Center of mass. + @param TW Tanks weights. + @param VOLS List of tank volumes. @param tanks Considered tanks. @param roll Roll angle. @param var_trim True if the trim angle should be recomputed at each roll angle, False otherwise. - @return GZ value + @return GZ value, equilibrium draft, and equilibrium trim angle (0 if + variable trim has not been requested) """ gz = 0.0 @@ -107,17 +118,21 @@ def solve_point(W, COG, TW, ship, tanks, roll, var_trim=True): max_disp / 1000.0 / G, (W + TW) / 1000.0 / G)) return None trim = 0.0 + for i in range(MAX_EQUILIBRIUM_ITERS): # Get the displacement, and the bouyance application point disp, B, Cb = Hydrostatics.displacement(ship, draft, roll, trim) disp *= 1000.0 * G # Add the tanks effect on the center of gravity - # TODO - # --- - + cog = Vector(COG.x * W, COG.y * W, COG.z * W) + for i,t in enumerate(tanks): + tank_weight = VOLS[i] * t[1] * G + cog += t[0].Proxy.getCoG(t[0], VOLS[i], roll, trim).multiply( + tank_weight / Units.Metre.Value) + cog = cog.multiply(1.0 / (W + TW)) # Compute the errors draft_error = -(disp - W - TW) / max_disp - R = COG - B + R = cog - B if not var_trim: trim_error = 0.0 else: @@ -131,9 +146,7 @@ def solve_point(W, COG, TW, ship, tanks, roll, var_trim=True): draft += draft_error * max_draft trim += trim_error - App.Console.PrintMessage("draft and trim: {}, {}\n".format(draft, trim)) - # GZ should be provided in the Free surface oriented frame of reference c = math.cos(math.radians(roll)) s = math.sin(math.radians(roll)) - return c * R.y - s * R.z \ No newline at end of file + return c * R.y - s * R.z, draft, trim \ No newline at end of file