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/AUVSim/.gitignore b/AUVSim/.gitignore new file mode 100644 index 0000000..ba0430d --- /dev/null +++ b/AUVSim/.gitignore @@ -0,0 +1 @@ +__pycache__/ \ No newline at end of file diff --git a/AUVSim/__init__.py b/AUVSim/__init__.py new file mode 100644 index 0000000..5c439b7 --- /dev/null +++ b/AUVSim/__init__.py @@ -0,0 +1,5 @@ +from .auv import ( + AUV as AUV, + AUVPath as AUVPath, + AUVPathGen as AUVPathGen, +) \ No newline at end of file diff --git a/AUVSim/auv.py b/AUVSim/auv.py new file mode 100644 index 0000000..73c3bb1 --- /dev/null +++ b/AUVSim/auv.py @@ -0,0 +1,65 @@ +from Raytrace import TriangleMesh +from Raytrace.SideScan import SideScan +import numpy as np +from typing import Any +from CustomTypes import Ray, Vector3, Vector3Array + +class AUV: + ideal_distance:float # Currently unused + start_pos:Vector3 + start_facing:Vector3 + def __init__(self, + position:Vector3, + facing:Vector3, + ideal_distance:float = -1, + ): + self.start_pos = position + self.start_facing = facing + self.ideal_distance = ideal_distance +class AUVPath: + positions:Vector3Array + facings:Vector3Array + + def __init__(self, positions:Vector3Array, facings:Vector3Array): + self.positions = positions + self.facings = facings + @property + def rays(self) -> list[Ray]: + return list(map(Ray, self.facings, self.positions)) + +class AUVPathGen: + mesh:TriangleMesh + auv:AUV + def __init__(self, mesh:TriangleMesh, auv:AUV): + self.mesh = mesh + self.auv = auv + def get_path(self, travel_distance:float, samples:int) -> AUVPath: + travel_step = travel_distance / (samples - 1) + positions = [self.auv.start_pos] + facings = [self.auv.start_facing] + + # Utility functions + # HACK: There should be a function to generate one side ray + # instead of generating a list with length one. + get_side_ray = lambda face, pos:\ + SideScan.generate_rays(Ray(face, pos), 0, 1, 2)[0] + get_dist = lambda ray: self.mesh.raytrace(ray) + normalize = lambda vector: vector / np.sqrt(np.dot(vector, vector)) + + for _ in range(samples - 1): + # find distance at current state + prev_dist = get_dist(get_side_ray(facings[-1], positions[-1])) + # move frowards + positions.append(positions[-1] + facings[-1] * travel_step) + # shoot side ray + # NOTE: facing has not yet been updated + side_ray = get_side_ray(facings[-1], positions[-1]) + ray_dist = self.mesh.raytrace(side_ray) + # point in the direction that the hull is sloping + rise = (ray_dist - prev_dist) * side_ray.direction + if not (np.isfinite(ray_dist) and np.isfinite(prev_dist)): + rise = Vector3(0, 0, 0) + new_facing = normalize(facings[-1] + rise) + facings.append(new_facing) + + return AUVPath(Vector3Array(positions), Vector3Array(facings)) \ No newline at end of file diff --git a/CustomTypes/.gitignore b/CustomTypes/.gitignore new file mode 100644 index 0000000..ba0430d --- /dev/null +++ b/CustomTypes/.gitignore @@ -0,0 +1 @@ +__pycache__/ \ No newline at end of file diff --git a/CustomTypes/__init__.py b/CustomTypes/__init__.py new file mode 100644 index 0000000..ff9b0ea --- /dev/null +++ b/CustomTypes/__init__.py @@ -0,0 +1,11 @@ +from .custom_np import CustomNP as CustomNP +from .custom_types import ( + Vector3 as Vector3, + Triangle as Triangle, + Ray as Ray, +) +from .custom_arrays import ( + FloatArray as FloatArray, + TriangleArray as TriangleArray, + Vector3Array as Vector3Array, +) \ No newline at end of file diff --git a/CustomTypes/custom_arrays.py b/CustomTypes/custom_arrays.py new file mode 100644 index 0000000..e881f0b --- /dev/null +++ b/CustomTypes/custom_arrays.py @@ -0,0 +1,54 @@ +from typing import Any, Iterable, Literal, overload +import numpy as np + +from . import CustomNP, Triangle, Vector3 + +class FloatArray(CustomNP): + def __new__(cls, + content: Iterable[float]|np.ndarray|None = None, + ): + if isinstance(content, np.ndarray): + if content.ndim != 1: raise ValueError(content.shape) # must be 1d + return content.view(cls) + if content is None: + return np.ndarray((0,), float).view(cls) + return np.array(content).view(cls) + def __getitem__(self, index): + if isinstance(index, slice): return super().__getitem__(index).view(type(self)) + return super().__getitem__(index) + +class Vector3Array(CustomNP): + def __new__(cls, + content: Iterable[Vector3]|np.ndarray|None = None, + ): + if isinstance(content, np.ndarray): + if content.ndim == 1: content = content.reshape((1,) + content.shape) + if content.ndim != 2: raise ValueError(content.shape) # must be 1d or 2d + if content.shape[1] != 3: raise ValueError(content.shape) # must be shape (n, 3) + return content.view(cls) + if content is None: + return np.ndarray((0, 3), float).view(cls) + return np.array(content).view(cls) + def __getitem__(self, index): + if isinstance(index, slice): return super().__getitem__(index).view(type(self)) + return super().__getitem__(index) + +class TriangleArray(CustomNP): + def __new__(cls, + content: Iterable[Triangle]|np.ndarray|None = None, + ): + if isinstance(content, np.ndarray): + if content.ndim == 2: content = content.reshape((1,) + content.shape) + if content.ndim != 3: raise ValueError(content.shape) # must be 2d or 3d + if content.shape[1] != 3 or content.shape[2] != 3: raise ValueError(content.shape) # must be shape (n, 3, 3) + return content.view(cls) + if content is None: + return np.ndarray((0, 3, 3), float).view(cls) + return np.array(content).view(cls) + def __getitem__(self, index): + if isinstance(index, slice): return super().__getitem__(index).view(type(self)) + return super().__getitem__(index) + @property + def verticies(self) -> tuple[Vector3Array, Vector3Array, Vector3Array]: + v0, v1, v2 = self.swapaxes(0, 1) + return (v0.view(Vector3Array), v1.view(Vector3Array), v2.view(Vector3Array)) \ No newline at end of file diff --git a/CustomTypes/custom_arrays.pyi b/CustomTypes/custom_arrays.pyi new file mode 100644 index 0000000..2b29ba9 --- /dev/null +++ b/CustomTypes/custom_arrays.pyi @@ -0,0 +1,53 @@ +import numpy as np +from typing import Any, Iterable, Literal, Self, Type, TypeVar, overload + +from . import CustomNP, Triangle, Vector3 + +class FloatArray(CustomNP[tuple[int, Literal[3]], np.dtype[np.float_]]): + @overload + def __new__(cls) -> FloatArray:... + @overload + def __new__(cls, content: Iterable[float]) -> FloatArray:... + @overload + def __new__(cls, content: np.ndarray[tuple[int], np.dtype[np.float_]]) -> FloatArray:... + + @overload # type: ignore + def __getitem__(self, index:int) -> float: ... + @overload + def __getitem__(self, index:slice) -> Self: ... + +class Vector3Array(CustomNP[tuple[int, Literal[3]], np.dtype[np.float_]]): + @overload + def __new__(cls) -> Vector3Array:... + @overload + def __new__(cls, content: Iterable[Vector3]) -> Vector3Array:... + @overload + def __new__(cls, content: np.ndarray[tuple[int, Literal[3]], np.dtype[np.float_]]) -> Vector3Array:... + + @overload # type: ignore + def __getitem__(self, index:int) -> Vector3: ... + @overload + def __getitem__(self, index:tuple[int, Literal[0, 1, 2]]) -> float: ... + @overload + def __getitem__(self, index:slice) -> Self: ... + +class TriangleArray(CustomNP[tuple[int, Literal[3], Literal[3]], np.dtype[np.float_]]): + @overload + def __new__(cls) -> TriangleArray:... + @overload + def __new__(cls, content: Iterable[Triangle]) -> TriangleArray:... + @overload + def __new__(cls, content: np.ndarray[tuple[int, Literal[3], Literal[3]], np.dtype[np.float_]]) -> TriangleArray:... + + @overload # type: ignore + def __getitem__(self, index:int) -> Triangle: ... + @overload + def __getitem__(self, index:tuple[int, Literal[0, 1, 2]]) -> Vector3: ... + @overload + def __getitem__(self, index:tuple[int, Literal[0, 1, 2], Literal[0, 1, 2]]) -> float: ... + @overload + def __getitem__(self, index:slice) -> Self: ... + + @property + def verticies(self) -> tuple[Vector3Array, Vector3Array, Vector3Array]: ... + diff --git a/CustomTypes/custom_np.py b/CustomTypes/custom_np.py new file mode 100644 index 0000000..9ab0a44 --- /dev/null +++ b/CustomTypes/custom_np.py @@ -0,0 +1,21 @@ + +from typing import Generic, TypeVar, Self +import numpy as np + + +_S = TypeVar('_S', bound=tuple) +_T = TypeVar('_T', bound=np.dtype) +class CustomNP(np.ndarray[_S, _T]): + def __add__(self, other): + return super().__add__(other).view(type(self)) + def __sub__(self, other): + return super().__sub__(other).view(type(self)) + def __mul__(self, other): + return super().__mul__(other).view(type(self)) + def __truediv__(self, other): + return super().__truediv__(other).view(type(self)) + def __rtruediv__(self, other): + return super().__rtruediv__(other).view(type(self)) + @property + def nd(self): + return self.view(np.ndarray) \ No newline at end of file diff --git a/CustomTypes/custom_np.pyi b/CustomTypes/custom_np.pyi new file mode 100644 index 0000000..9786798 --- /dev/null +++ b/CustomTypes/custom_np.pyi @@ -0,0 +1,15 @@ +from typing import TypeVar, Self +import numpy as np + + +_S = TypeVar('_S', bound=tuple) +_T = TypeVar('_T', bound=np.dtype) +class CustomNP(np.ndarray[_S, _T]): + def __add__(self, other:Self) -> Self: ... # type: ignore + def __sub__(self, other:Self) -> Self: ... # type: ignore + def __mul__(self, other:Self|float) -> Self: ... # type: ignore + def __truediv__(self, other:Self|float) -> Self: ... # type: ignore + def __rtruediv__(self, other:float) -> Self: ... # type: ignore + + @property + def nd(self) -> np.ndarray: ... diff --git a/CustomTypes/custom_types.py b/CustomTypes/custom_types.py new file mode 100644 index 0000000..0057c89 --- /dev/null +++ b/CustomTypes/custom_types.py @@ -0,0 +1,42 @@ +from typing import Any, Iterable, Literal, Self, overload +import numpy as np +from . import CustomNP + +class Vector3(CustomNP): + def __new__(cls, + x: float | tuple[float, float, float] | Self, y = 0, z = 0): + if isinstance(x, tuple): + x, y, z = x + elif isinstance(x, Vector3): + x, y, z = x + assert isinstance(x, (float, int, np.floating)) + return np.array([x, y, z]).view(cls) + def __str__(self): + return f"({self[0]}, {self[1]}, {self[2]})" + def __repr__(self): + return f"Vector3({self[0]}, {self[1]}, {self[2]})" + +class Triangle(CustomNP): + def __new__(cls, + v0: Vector3 | tuple[Vector3, Vector3, Vector3], + v1: Vector3 = Vector3(0, 0, 0), + v2: Vector3 = Vector3(0, 0, 0), + ) -> 'Triangle': + if isinstance(v0, tuple): + v0, v1, v2 = v0 + return np.array([v0, v1, v2]).view(cls) + + def __str__(self) -> str: + return f"({self[0]}, {self[1]}, {self[2]})" + def __repr__(self) -> str: + return f"Vector3({self[0]}, {self[1]}, {self[2]})" + +class Ray: + direction:Vector3 + origin:Vector3 + def __init__(self, + direction:Vector3, + origin:Vector3, + ) -> None: + self.direction = direction + self.origin = origin \ No newline at end of file diff --git a/CustomTypes/custom_types.pyi b/CustomTypes/custom_types.pyi new file mode 100644 index 0000000..04fac9a --- /dev/null +++ b/CustomTypes/custom_types.pyi @@ -0,0 +1,40 @@ +from typing import Any, Iterable, Iterator, Literal, TypeAlias, TypeVar, overload +import numpy as np +from . import CustomNP + +class Vector3(CustomNP[tuple[Literal[3]], np.dtype[np.float_]]): + @overload + def __new__(cls, x: float, y: float, z: float) -> Vector3: ... + @overload + def __new__(cls, x: tuple[float, float, float]) -> Vector3: ... + @overload + def __new__(cls, x: Vector3) -> Vector3: ... + + def __getitem__(self, index:Literal[0, 1, 2]) -> float: ... # type: ignore + +class Triangle(CustomNP[tuple[Literal[3], Literal[3]], np.dtype[np.float_]], Iterable[Vector3]): + @overload + def __new__(cls, v0: Vector3, v1: Vector3, v2: Vector3) -> Triangle: ... + @overload + def __new__(cls, v0: tuple[Vector3, Vector3, Vector3]) -> Triangle: ... + @overload + def __new__(cls, v0: Triangle) -> Triangle: ... + + @overload # type: ignore + def __getitem__(self, index:Literal[0, 1, 2]) -> Vector3: ... + @overload + def __getitem__(self, index:tuple[Literal[0, 1, 2], Literal[0, 1, 2]]) -> float: ... + + def __add__(self, other:Triangle) -> Triangle: ... # type: ignore + def __sub__(self, other:Triangle) -> Triangle: ... # type: ignore + + def __iter__(self) -> Iterator[Vector3]: ... +class Ray: + direction:Vector3 + origin:Vector3 + def __init__(self, + direction:Vector3, + origin:Vector3, + ) -> None: + self.direction = direction + self.origin = origin \ No newline at end of file diff --git a/Generate/.gitignore b/Generate/.gitignore new file mode 100644 index 0000000..ba0430d --- /dev/null +++ b/Generate/.gitignore @@ -0,0 +1 @@ +__pycache__/ \ No newline at end of file diff --git a/Generate/__init__.py b/Generate/__init__.py new file mode 100644 index 0000000..94cac17 --- /dev/null +++ b/Generate/__init__.py @@ -0,0 +1,13 @@ +from .auv import( + place_auvs_over_points as place_auvs_over_points, + generate_paths as generate_paths, + scan_paths as scan_paths, +) +from .hull import ( + generate_hulls as generate_hulls, + generate_points_on_hull as generate_points_on_hull, + generate_vertical_bounds as generate_vertical_bounds, + load_anomalies as load_anomalies, + pick_anomalies as pick_anomalies, + place_anomalies_at_points as place_anomalies_at_points, +) \ No newline at end of file diff --git a/Generate/auv.py b/Generate/auv.py new file mode 100644 index 0000000..9a82f86 --- /dev/null +++ b/Generate/auv.py @@ -0,0 +1,114 @@ +from typing import Callable + +from CustomTypes import Vector3, Vector3Array +from Raytrace import TriangleMesh +from AUVSim import AUV, AUVPathGen, AUVPath +from Raytrace.SideScan import SideScan, ScanReading, RawScanReading + +from .helper import get_surface_normal, normalize + +def place_auv_over_point( + point:Vector3, + mesh:TriangleMesh, + distance: float, + align_start_facing = True, + use_cw_normal:None|bool = None, + ) -> AUV: + normal = get_surface_normal(mesh, point, use_cw_normal) + position = point + normal * distance + facing = Vector3(1, 0, 0) + if align_start_facing: + facing = normalize(Vector3(-normal[1], normal[0], 0)) # tangent w/ no z + + return AUV(position, facing, distance) + +def place_auvs_over_points( + points:Vector3Array, + meshs:TriangleMesh|list[TriangleMesh], + distance: float, + align_start_facing = True, + use_cw_normal:None|bool = None, + ) -> list[AUV]: + if not isinstance(meshs, list): + meshs = [meshs] * len(points) + return [ + place_auv_over_point(point, mesh, distance, align_start_facing, use_cw_normal) + for point, mesh in zip(points, meshs) + ] + +def generate_path( + mesh:TriangleMesh, + auv:AUV, + travel_distance:float, + samples:int, + back_steps:int|None = None, + ) -> AUVPath: + # Step 1: Backstep + if not (back_steps is None): + auv_copy = AUV(auv.start_pos, auv.start_facing) + auv_copy.ideal_distance = auv.ideal_distance + back_dist = travel_distance / (samples - 1) * back_steps + auv = step_back_auv(mesh, auv_copy, back_dist, back_steps) + # Step 2: Get path + return AUVPathGen(mesh, auv).get_path(travel_distance, samples) + +def generate_paths( + meshs:TriangleMesh|list[TriangleMesh], + auvs:list[AUV], + travel_distance:float, + samples:int, + back_steps:int|None|list[int|None] = None, + ) -> list[AUVPath]: + if not isinstance(meshs, list): + meshs = [meshs] * len(auvs) + if not isinstance(back_steps, list): + back_steps = [back_steps] * len(auvs) + return [ + generate_path(mesh, auv, travel_distance, samples, bs) + for mesh, auv, bs in zip(meshs, auvs, back_steps) + ] + +def step_back_auv( + mesh:TriangleMesh, + auv:AUV, + travel_distance:float, + steps:int, + ) -> AUV: + auv.start_facing = auv.start_facing * -1 + path_gen = AUVPathGen(mesh, auv) + back_start = path_gen.get_path(travel_distance, steps + 1).rays[-1] + auv.start_facing = auv.start_facing * -1 # undo changes + + return AUV(back_start.origin, back_start.direction * -1, auv.ideal_distance) + +def scan_path( + scanner:SideScan, + path: AUVPath, + min_angle:float, + max_angle:float, + angle_reselution:int, + silent:bool = False, + process:Callable[[RawScanReading], ScanReading] = ScanReading + ) -> list[ScanReading]: + readings:list[ScanReading] = [] + for n, orientation in enumerate(path.rays): + rays = SideScan.generate_rays(orientation, min_angle, max_angle, angle_reselution) + if not silent: + print(n + 1, len(path.rays), sep='/', end=' \r') + readings.append(process(scanner.scan_rays(rays))) + return readings +def scan_paths( + scanners:SideScan|list[SideScan], + paths: list[AUVPath], + min_angle:float, + max_angle:float, + angle_reselution:int, + silent:bool = False, + process:Callable[[RawScanReading], ScanReading] = ScanReading + ) -> list[list[ScanReading]]: + if not isinstance(scanners, list): + scanners = [scanners] * len(paths) + return [ + scan_path(scanner, path, min_angle, max_angle, angle_reselution, silent, process) + for scanner, path in zip(scanners, paths) + ] diff --git a/Generate/helper.py b/Generate/helper.py new file mode 100644 index 0000000..d255759 --- /dev/null +++ b/Generate/helper.py @@ -0,0 +1,38 @@ +import numpy as np + +from CustomTypes import Ray, Vector3 +from Raytrace import TriangleMesh + + +def normalize(v:Vector3) -> Vector3: + return v / float(np.sqrt(v.dot(v))) + +def get_surface_normal( + mesh:TriangleMesh, + point:Vector3, + use_cw_normal:None|bool = None, + ): + # Step 1: Find the triangle at the point + # NOTE: Using the closest centroid is not necesarily correct + # But its a decent aproximation in most cases + min_ind = find_closest_centroid_ind(mesh, point) + tri = mesh.triangles[min_ind] + if tri is None: raise ValueError(f"sample_ray could not find surface at point {point}") + # Step 2: Find the normal at the triangle + normal = np.cross(tri[1] - tri[0], tri[2] - tri[0]).view(Vector3) + normal = normalize(normal) + if use_cw_normal: + normal = normal * -1 + if use_cw_normal is None: # Guess at correct normal direction + if (normal[1] < 0) != (point[1] < 0): + normal = normal * -1 # Point normal away from xz plane + return normal + +def find_closest_centroid_ind( + mesh:TriangleMesh, + point:Vector3, + ) -> int: + # NOTE: Could (probably) be sped up with scipy.spatial.KDTree + rel_pos = mesh.centroids - point + sqr_dist = np.sum(rel_pos*rel_pos, axis=-1) + return int(np.argmin(sqr_dist)) \ No newline at end of file diff --git a/Generate/hull.py b/Generate/hull.py new file mode 100644 index 0000000..55a33aa --- /dev/null +++ b/Generate/hull.py @@ -0,0 +1,122 @@ +from copy import deepcopy +from typing import overload + +import numpy as np +from CustomTypes import Vector3, TriangleArray, Vector3Array +from Raytrace import TriangleMesh, CompositeMesh, BVHMesh +import ShipGen as SG +import random + +def generate_hulls( + find_best:int, + out_of:int, + x_length:float|None = None, + characteristics:SG.ShipCharacteristics|None = None, + silent:bool = False, + ) -> list[BVHMesh]: + # Step 1: Generate the hulls + raw_hulls = SG.generate_hulls(find_best, out_of, characteristics) + # Step 2: Apply scale + if not (x_length is None): + for hull in raw_hulls: + SG.normalize_x(hull, x_length) + # Step 3: Convert them to BVHs + triangles_list = SG.mesh_list_to_triangles(raw_hulls) + hulls:list[BVHMesh] = [] + for n, triangles in enumerate(triangles_list): + if not silent: + print(f'Building BVH {n} / {len(triangles_list)}') + hulls.append(BVHMesh(triangles, min_node_size = 100)) + return hulls + +def place_anomaly_at_point( + hull:TriangleMesh, + anomaly:TriangleMesh, + point:Vector3, + ) -> CompositeMesh: + anomaly = deepcopy(anomaly) + # TODO: A proper addition method should be implemented + # HACK: Will only work fully for regular TriangleMesh + if type(anomaly) != TriangleMesh: + raise NotImplementedError('place_anomaly_at_point only works for anomalies of type TriangleMesh (no subclasses)') + anomaly.triangles = TriangleArray(anomaly.triangles + point) + + return CompositeMesh([hull, anomaly]) +def place_anomalies_at_points( + hull:TriangleMesh, + anomalies:list[TriangleMesh], + points:Vector3Array, + ) -> list[CompositeMesh]: + return [ + place_anomaly_at_point(hull, anomaly, point) + for anomaly, point in zip(anomalies, points) + ] + +def generate_point_on_hull_anywhere( + hull:TriangleMesh, + ) -> Vector3: + # NOTE: not a perfectly even distribution + triangle_ind = random.randint(0, len(hull.triangles)) + triangle = hull.triangles[triangle_ind] + rand_x = random.random() + rand_y = random.random() + if rand_y > 1 - rand_x: # outside triangle + rand_x = 1 - rand_x + rand_y = 1 - rand_y + off_x = (triangle[1] - triangle[0]) * rand_x + off_y = (triangle[2] - triangle[0]) * rand_y + return triangle[0] + off_x + off_y +def generate_point_on_hull_in( + hull:TriangleMesh, + min_bound:Vector3, + max_bound:Vector3, + attempts:int = 1000 + ) -> Vector3: + for _ in range(attempts): + point = generate_point_on_hull_anywhere(hull) + if np.all(min_bound <= point) and np.all(point <= max_bound): + return point + raise RuntimeError(f'Was not able to generate valid point within {attempts} attempts') + +def generate_point_on_hull( + hull:TriangleMesh, + min_bound:Vector3|None = None, + max_bound:Vector3|None = None, + attempts:int = 1000 + ) -> Vector3: + if min_bound is None or max_bound is None: + return generate_point_on_hull_anywhere(hull) + else: + return generate_point_on_hull_in(hull, min_bound, max_bound, attempts) +def generate_points_on_hull( + count:int, + hull:TriangleMesh, + min_bound:Vector3|None = None, + max_bound:Vector3|None = None, + attempts:int = 1000 + ) -> Vector3Array: + return Vector3Array([ + generate_point_on_hull(hull, min_bound, max_bound, attempts) + for _ in range(count)]) + +def generate_vertical_bounds( + hull:TriangleMesh, + cut_percentage:float, + ) -> tuple[Vector3, Vector3]: + min_z = np.min(hull.triangles.nd[:,2]) + max_z = np.max(hull.triangles.nd[:,2]) + height = max_z - min_z + min_bound = Vector3(-np.inf, -np.inf, min_z + height * cut_percentage) + max_bound = Vector3( np.inf, np.inf, max_z - height * cut_percentage) + return min_bound, max_bound + +def load_anomalies(paths:str|list[str]) -> list[TriangleMesh]: + if isinstance(paths, str): return load_anomalies([paths]) + return [TriangleMesh(stl_path = path) for path in paths] +def pick_anomalies( + count:int, + anomalies:list[TriangleMesh], + weights:list[float], + ) -> list[TriangleMesh]: + #TODO: Implement weights + return [random.choice(anomalies) for _ in range(count)] \ No newline at end of file diff --git a/Generate/outline.ipynb b/Generate/outline.ipynb new file mode 100644 index 0000000..a5933cd --- /dev/null +++ b/Generate/outline.ipynb @@ -0,0 +1,167 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# This notebook walks through how to create items for the dataset" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from typing import Callable\n", + "import sys; sys.path.append('..')\n", + "\n", + "# from Generate import (\n", + "# generate_hulls,\n", + "# generate_paths,\n", + "# place_auvs_over_points,\n", + "# generate_points_on_hull,\n", + "# generate_vertical_bounds,\n", + "# load_anomalies,\n", + "# pick_anomalies,\n", + "# place_anomalies_at_points,\n", + "# scan_paths,\n", + "# )\n", + "from Generate import *\n", + "\n", + "from ShipGen import ShipCharacteristics\n", + "from Raytrace.SideScan import SideScan, RawScanReading, ScanReading" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "### Parameters ###\n", + "hull_count = 10\n", + "reject_hulls = 40\n", + "reg_readings_per_hull = 2\n", + "anom_readings_per_hull = 2\n", + "anomalies = load_anomalies() # TODO\n", + "ideal_distance = 1\n", + "scan_distance = 1\n", + "readings_per_scan = 100\n", + "rays_per_reading = 1000\n", + "result_reselution = 1000\n", + "min_angle = np.pi/4\n", + "max_angle = -np.pi/2\n", + "\n", + "randomize_reading:Callable[[RawScanReading], ScanReading] = ScanReading\n", + "\n", + "# Less imortant\n", + "x_length = 10\n", + "hull_characteristics:ShipCharacteristics|None = None\n", + "vertical_deadzone_percent = 0.1\n", + "smooth_dist = 0.05" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "### Generate hulls ###\n", + "hulls = generate_hulls(hull_count, hull_count + reject_hulls, x_length, hull_characteristics)\n", + "# Loop over the hulls (only showing one for the sake of example)\n", + "# for hull in hulls:\n", + "hull = hulls[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "### Find points on hulls ###\n", + "min_b, max_b = generate_vertical_bounds(hull, vertical_deadzone_percent)\n", + "total_points = reg_readings_per_hull + anom_readings_per_hull\n", + "points_combined = generate_points_on_hull(total_points, hull, min_b, max_b)\n", + "reg_points = points_combined[:reg_readings_per_hull]\n", + "anom_points = points_combined[reg_readings_per_hull:]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "### Place anomalies ###\n", + "anomaly_order = pick_anomalies(anom_readings_per_hull, anomalies, [])\n", + "anomaly_hulls = place_anomalies_at_points(hull, anomaly_order, anom_points)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "### Place auvs ###\n", + "reg_auvs = place_auvs_over_points(reg_points, hull, ideal_distance)\n", + "anom_auvs = place_auvs_over_points(anom_points, anomaly_hulls, ideal_distance)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "### Scan hulls w/wo anomalies ###\n", + "# TODO: add ranomization to backsteps\n", + "back_steps = readings_per_scan // 2\n", + "reg_paths = generate_paths(hull, reg_auvs, scan_distance, readings_per_scan, back_steps)\n", + "scanner = SideScan(hull, smooth_dist, result_reselution)\n", + "reg_scans = scan_paths(scanner, reg_paths, min_angle, max_angle, rays_per_reading, process = randomize_reading)\n", + "\n", + "anom_paths = generate_paths(anomaly_hulls, anom_auvs, scan_distance, readings_per_scan, back_steps)\n", + "anom_scanners = [SideScan(h, smooth_dist, result_reselution) for h in anomaly_hulls]\n", + "anom_scans = scan_paths(anom_scanners, anom_paths, min_angle, max_angle, rays_per_reading, process = randomize_reading)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "### Done ###\n", + "reg_scans, anom_scans" + ] + } + ], + "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/Raytrace/CompositeMesh.py b/Raytrace/CompositeMesh.py deleted file mode 100644 index f1dc5b6..0000000 --- a/Raytrace/CompositeMesh.py +++ /dev/null @@ -1,11 +0,0 @@ -import numpy as np -from Raytrace.TriangleMesh import TriangleMesh, Ray - -class CompositeMesh(TriangleMesh): - meshes: list[TriangleMesh] - def __init__(self, meshes:list[TriangleMesh]) -> None: - self.meshes = meshes - self.triangles = np.concatenate([mesh.triangles for mesh in self.meshes]) - def raytrace(self, ray:Ray) -> float: - distances = [mesh.raytrace(ray) for mesh in self.meshes] - return min(distances, default = np.inf) \ No newline at end of file diff --git a/Raytrace/Examples/ExampleMultipleSonarReadings.ipynb b/Raytrace/Examples/ExampleMultipleSonarReadings.ipynb index d5f6cb9..9d5f4e9 100644 --- a/Raytrace/Examples/ExampleMultipleSonarReadings.ipynb +++ b/Raytrace/Examples/ExampleMultipleSonarReadings.ipynb @@ -2,13 +2,13 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sys; sys.path.append('../..')\n", - "from Raytrace.TriangleMesh import Ray\n", - "from Raytrace.BVHMesh import BVHMesh\n", + "from CustomTypes import Ray, Vector3, Vector3Array\n", + "from Raytrace import BVHMesh\n", "from Raytrace.SideScan import SideScan, ScanReading\n", "import numpy as np\n", "import time" @@ -34,18 +34,18 @@ "metadata": {}, "outputs": [], "source": [ - "facing = np.array([1, 0, 0])\n", + "facing = Vector3(1, 0, 0)\n", "min_angle = 0\n", "max_angle = -np.pi/2\n", "sample_ray_count = 1000\n", "\n", - "origin_start = np.array([5, -1, 1])\n", - "origin_end = np.array([9, -1, 1])\n", + "origin_start = Vector3(5, -1, 1)\n", + "origin_end = Vector3(9, -1, 1)\n", "readings_count = 100\n", "# liniar interpolation between the two for this simple demo\n", - "origins = [origin_start * (1-i) + origin_end * i for i in np.arange(0, 1+1/(readings_count-1)/2, 1/(readings_count-1))]\n", + "origins = Vector3Array([origin_start * (1-i) + origin_end * i for i in np.arange(0, 1+1/(readings_count-1)/2, 1/(readings_count-1))])\n", "\n", - "orientations = [Ray(facing, origin) for origin in origins]\n", + "orientations = list(map(Ray, [facing] * readings_count, origins))\n", "\n", "rays_list = [SideScan.generate_rays(orientation, min_angle, max_angle, sample_ray_count) for orientation in orientations]\n", "\n", diff --git a/Raytrace/Examples/ExampleSingleSonarReading.ipynb b/Raytrace/Examples/ExampleSingleSonarReading.ipynb index 0a89788..4c4bed6 100644 --- a/Raytrace/Examples/ExampleSingleSonarReading.ipynb +++ b/Raytrace/Examples/ExampleSingleSonarReading.ipynb @@ -2,16 +2,17 @@ "cells": [ { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sys; sys.path.append('../..')\n", - "from Raytrace.TriangleMesh import Ray\n", - "from Raytrace.BVHMesh import BVHMesh\n", + "from CustomTypes import Ray\n", + "from Raytrace import BVHMesh\n", "from Raytrace.SideScan import SideScan, ScanReading\n", "import numpy as np\n", - "import time" + "import time\n", + "from CustomTypes import Vector3" ] }, { @@ -34,8 +35,8 @@ "metadata": {}, "outputs": [], "source": [ - "origin = np.array([5, -1, 1])\n", - "facing = np.array([1, 0, 0])\n", + "origin = Vector3(5, -1, 1)\n", + "facing = Vector3(1, 0, 0)\n", "min_angle = 0\n", "max_angle = -np.pi/2\n", "sample_ray_count = 1000\n", @@ -47,7 +48,7 @@ "raw_reading = SideScan(mesh).scan_rays(rays)\n", "reading = ScanReading(raw_reading)\n", "\n", - "print('Triangles:', mesh.triangles.shape[0])\n", + "print('Triangles:', len(mesh.triangles))\n", "raw_reading.print_summary()" ] }, diff --git a/Raytrace/Meshs/__init__.py b/Raytrace/Meshs/__init__.py new file mode 100644 index 0000000..cc039d9 --- /dev/null +++ b/Raytrace/Meshs/__init__.py @@ -0,0 +1,3 @@ +from .triangle_mesh import TriangleMesh as TriangleMesh +from .bvh_mesh import BVHMesh as BVHMesh +from .composite_mesh import CompositeMesh as CompositeMesh \ No newline at end of file diff --git a/Raytrace/BVHMesh.py b/Raytrace/Meshs/bvh_mesh.py similarity index 71% rename from Raytrace/BVHMesh.py rename to Raytrace/Meshs/bvh_mesh.py index 6edb52e..14dce36 100644 --- a/Raytrace/BVHMesh.py +++ b/Raytrace/Meshs/bvh_mesh.py @@ -1,34 +1,37 @@ -from typing import Self, Literal +from typing import Self, Literal, Type, TypeVar import numpy as np -from Raytrace.TriangleMesh import TriangleMesh, Ray + +from . import TriangleMesh +from CustomTypes import Ray, TriangleArray, Vector3, Vector3Array, Triangle # Adapted from https://jacco.ompf2.com/2022/04/13/how-to-build-a-bvh-part-1-basics/ class BVHAABB: - min_pos:np.ndarray - max_pos:np.ndarray + min_pos:Vector3 + max_pos:Vector3 def __init__(self): - self.min_pos = np.array([np.nan]*3) - self.max_pos = np.array([np.nan]*3) + self.min_pos = Vector3((np.nan,)*3) + self.max_pos = Vector3((np.nan,)*3) @staticmethod - def from_array(triangles:np.ndarray) -> 'BVHAABB': + def from_array(triangles:TriangleArray) -> 'BVHAABB': out = BVHAABB() - if triangles.size == 0: return out - triangles = triangles.reshape((-1, 3)) - out.min_pos = triangles.min(axis = 0) - out.max_pos = triangles.max(axis = 0) + if len(triangles) == 0: return out + verts = triangles.reshape((-1, 3)).view(Vector3Array) + out.min_pos = verts.min(axis = 0).view(Vector3) + out.max_pos = verts.max(axis = 0).view(Vector3) + return out def raytrace(self, ray:Ray) -> float: old_error_state = np.seterr(divide='ignore') - t1 = (self.min_pos - ray.origin) * ray.inverse_direction - t2 = (self.max_pos - ray.origin) * ray.inverse_direction + t1 = (self.min_pos - ray.origin) / ray.direction + t2 = (self.max_pos - ray.origin) / ray.direction np.seterr(**old_error_state) t = np.stack([t1, t2]) tmin:float = t.min(axis=0).max() tmax:float = t.max(axis=0).min() return tmin if tmax >= tmin and tmax > 0 else np.inf - def grow(self, point:np.ndarray): - self.min_pos = np.nanmin(np.stack([point, self.min_pos]), axis=0) - self.max_pos = np.nanmax(np.stack([point, self.max_pos]), axis=0) + def grow(self, point:Vector3): + self.min_pos = np.nanmin(np.array([point, self.min_pos]), axis=0).view(Vector3) + self.max_pos = np.nanmax(np.array([point, self.max_pos]), axis=0).view(Vector3) def grow_aabb(self, other:Self): self.grow(other.min_pos) self.grow(other.max_pos) @@ -43,27 +46,25 @@ class BVHNode: right:Self start_index:int tri_count:int - def get_subarray(self, triangles:np.ndarray) -> np.ndarray: + _T = TypeVar('_T', TriangleArray, Vector3Array) + def get_subarray(self, triangles:_T) -> _T: return triangles[self.start_index:self.start_index + self.tri_count] - def update_bounds(self, triangles:np.ndarray): + def update_bounds(self, triangles:TriangleArray): self.aabb = BVHAABB.from_array(self.get_subarray(triangles)) def is_leaf(self) -> bool: return self.tri_count > 0 class BVHMesh(TriangleMesh): - centroids:np.ndarray root:BVHNode node_count:int def __init__(self, *args, min_node_size = 2, **kwargs) -> None: super().__init__(*args, **kwargs) self.build_BVH(min_node_size) def build_BVH(self, min_node_size = 100) -> None: - # calculate triangle centroids for partitioning - self.centroids = self.triangles.sum(axis = 1) / 3 # assign all triangles to root node self.node_count = 1 self.root = root = BVHNode() root.start_index = 0 - root.tri_count = self.triangles.shape[0] + root.tri_count = len(self.triangles) if root.tri_count == 0: return root.update_bounds(self.triangles) # subdivide recursively @@ -80,8 +81,9 @@ def subdivide(self, node:BVHNode, min_node_size = 100) -> None: while i <= j: if self.centroids[i, axis] < split_pos: i += 1 else: - self.triangles[[i,j]] = self.triangles[[j,i]] - self.centroids[[i,j]] = self.centroids[[j,i]] + self.triangles.nd[[i,j]] = self.triangles.nd[[j,i]] + self.centroids.nd[[i,j]] = self.centroids.nd[[j,i]] + j -= 1 # abort split if one of the sides is empty left_count = i - node.start_index @@ -103,18 +105,22 @@ def subdivide(self, node:BVHNode, min_node_size = 100) -> None: self.subdivide(left) self.subdivide(right) - def raytrace(self, ray:Ray) -> float: + def raytrace_as_tri(self, ray:Ray) -> tuple[float, Triangle|None]: return self.BVH_raytrace(ray) - def BVH_raytrace(self, ray:Ray) -> float: + def BVH_raytrace(self, ray:Ray) -> tuple[float, Triangle|None]: node:BVHNode = self.root stack:list[BVHNode] = [] - best:float = np.inf + best_dist:float = np.inf + best_tri:Triangle|None = None while 1: if node.is_leaf(): intersections = self.batch_triangle_ray_intersection(node.get_subarray(self.triangles), ray) intersections[intersections < 0] = np.inf - best = min(best, np.min(intersections)) + min_index = int(np.argmin(intersections)) + if intersections[min_index] < best_dist: + best_dist = intersections[min_index] + best_tri = node.get_subarray(self.triangles)[min_index] if len(stack) == 0: break node = stack.pop() continue @@ -131,25 +137,27 @@ def BVH_raytrace(self, ray:Ray) -> float: node = child1 if np.isfinite(dist2): stack.append(child2) - return best - def find_best_split(self, node:BVHNode) -> tuple[float, int, float]: + return best_dist, best_tri + + def find_best_split(self, node:BVHNode) -> tuple[float, Literal[0, 1, 2], float]: BINS = 8 - axis:int = -1 + axis:Literal[0, 1, 2] = 0 split_pos:float = 0 best_cost:float = np.inf - triangles:np.ndarray[tuple[int, Literal[3], Literal[3]], np.dtype[np.float32]] = node.get_subarray(self.triangles) - centroids:np.ndarray[tuple[int, Literal[3]], np.dtype[np.float32]] = node.get_subarray(self.centroids) - for axis_i in range(3): - bounds_min = float(np.min(centroids[:,axis_i])) - bounds_max = float(np.max(centroids[:,axis_i])) + triangles:TriangleArray = node.get_subarray(self.triangles) + centroids:Vector3Array = node.get_subarray(self.centroids) + for axis_i in (0, 1, 2): + bounds_min = float(np.min(centroids.nd[:,axis_i])) + bounds_max = float(np.max(centroids.nd[:,axis_i])) if bounds_min == bounds_max: continue # populate the bins scale:float = BINS / (bounds_max - bounds_min) - bin_idx = ((centroids[:,axis_i] - bounds_min) * scale).astype(int) + bin_idx = ((centroids.nd[:,axis_i] - bounds_min) * scale).astype(int) bin_idx[bin_idx > BINS - 1] = BINS - 1 bin_bounds = [BVHAABB.from_array(triangles[bin_idx == n]) for n in range(BINS)] bin_counts = [np.count_nonzero(bin_idx == n) for n in range(BINS)] + # gather data for the 7 planes between the 8 bins left_area = np.zeros(BINS - 1, float) right_area = np.zeros(BINS - 1, float) diff --git a/Raytrace/Meshs/composite_mesh.py b/Raytrace/Meshs/composite_mesh.py new file mode 100644 index 0000000..6292cb3 --- /dev/null +++ b/Raytrace/Meshs/composite_mesh.py @@ -0,0 +1,22 @@ +from functools import cached_property +import numpy as np +from . import TriangleMesh +from CustomTypes import Ray, TriangleArray, Triangle + +class CompositeMesh(TriangleMesh): + meshes: list[TriangleMesh] + + def __init__(self, meshes:list[TriangleMesh]) -> None: + self.meshes = meshes + + def raytrace_as_tri(self, ray:Ray) -> tuple[float, Triangle|None]: + traces = [mesh.raytrace_as_tri(ray) for mesh in self.meshes] + if len(traces) == 0: return (np.inf, None) + return min(traces, key= lambda v: v[0]) + + @property # only should be used for rendering, not calculations + def triangles(self) -> TriangleArray: + return np.concatenate([mesh.triangles for mesh in self.meshes]).view(TriangleArray) + @triangles.setter + def triangles(self, _) -> None: + raise AttributeError('CompositeMesh does not support writing to its triangles property') diff --git a/Raytrace/TriangleMesh.py b/Raytrace/Meshs/triangle_mesh.py similarity index 51% rename from Raytrace/TriangleMesh.py rename to Raytrace/Meshs/triangle_mesh.py index d1b8898..e8c31ab 100644 --- a/Raytrace/TriangleMesh.py +++ b/Raytrace/Meshs/triangle_mesh.py @@ -1,31 +1,30 @@ +from typing import Any import numpy as np import pickle - -class Ray: - direction:np.ndarray - origin:np.ndarray - inverse_direction:np.ndarray - def __init__(self, direction:np.ndarray, origin:np.ndarray) -> None: - self.direction = direction - self.origin = origin - old_error_state = np.seterr(divide='ignore', invalid='ignore') - self.inverse_direction = np.ones(3) / direction - np.seterr(**old_error_state) +from CustomTypes import Triangle, TriangleArray, Ray, Vector3, Vector3Array, FloatArray +import stl class TriangleMesh: - # An n by 3 by 3 array of floats representing a list of triangles in the form triangle_index, vertex_index, axis_index - triangles:np.ndarray - - def __init__(self, triangles: np.ndarray|None = None, stl_path: str|None = None): + triangles:TriangleArray + def __init__(self, + triangles: np.ndarray|TriangleArray|None = None, + stl_path: str|None = None): if not stl_path is None: - import stl mesh = stl.Mesh.from_file(stl_path) - triangles = np.array([mesh.v0, mesh.v1, mesh.v2]).swapaxes(0, 1) - self.triangles = np.ndarray((0,3), float) if triangles is None else triangles + self.triangles = TriangleArray(np.array([mesh.v0, mesh.v1, mesh.v2]).swapaxes(0, 1)) + elif triangles is None: + self.triangles = TriangleArray([]) + elif isinstance(triangles, np.ndarray): + self.triangles = TriangleArray(triangles) + else: + self.triangles = triangles def raytrace(self, ray:Ray) -> float: - return self.array_raytrace(ray) + return self.raytrace_as_tri(ray)[0] + def raytrace_as_tri(self, ray:Ray) -> tuple[float, Triangle|None]: + return self.array_raytrace(ray) + def simple_raytrace(self, ray:Ray) -> float: best_dist = np.inf @@ -38,40 +37,41 @@ def simple_raytrace(self, ray:Ray) -> float: if dist < best_dist: best_dist = dist return best_dist - def array_raytrace(self, ray:Ray) -> float: + def array_raytrace(self, ray:Ray) -> tuple[float, Triangle|None]: out = self.batch_triangle_ray_intersection(self.triangles, ray) out[out < 0] = np.inf - - return np.min(out) + min_index = int(np.argmin(out)) + if np.isinf(out[min_index]): return (np.inf, None) # hit nothing + return (out[min_index], self.triangles[int(np.argmin(out))]) @staticmethod - def triangle_ray_intersection(triangle, ray:Ray, epsilon = 1e-10) -> float: + def triangle_ray_intersection(triangle:Triangle, ray:Ray, epsilon = 1e-10) -> float: # Translated from https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm#C++_implementation v0, v1, v2 = triangle edge1 = v1 - v0 edge2 = v2 - v0 - ray_cross_e2 = np.cross(ray.direction, edge2) - det = np.dot(edge1, ray_cross_e2) + ray_cross_e2 = np.cross(ray.direction, edge2).view(Vector3) + det:float = np.dot(edge1, ray_cross_e2) if -epsilon < det < epsilon: return -1 # This ray is parallel to this triangle. inv_det = 1.0 / det s = ray.origin - v0 - u = inv_det * np.dot(s, ray_cross_e2) + u:float = inv_det * np.dot(s, ray_cross_e2) if u < -epsilon or 1 < u: return -1 - s_cross_e1 = np.cross(s, edge1) - v = inv_det * np.dot(ray.direction, s_cross_e1) + s_cross_e1 = np.cross(s, edge1).view(Vector3) + v:float = inv_det * np.dot(ray.direction, s_cross_e1) if v < -epsilon or 1 < u + v: return -1 # At this stage we can compute t to find out where the intersection point is on the line. - t = inv_det * np.dot(edge2, s_cross_e1) + t:float = inv_det * np.dot(edge2, s_cross_e1) if t > epsilon: # ray intersection return t @@ -79,28 +79,28 @@ def triangle_ray_intersection(triangle, ray:Ray, epsilon = 1e-10) -> float: else: # This means that there is a line intersection but not a ray intersection. return -1 @staticmethod - def batch_triangle_ray_intersection(triangle_array, ray:Ray, epsilon = 1e-10) -> np.ndarray: + def batch_triangle_ray_intersection(triangle_array:TriangleArray, ray:Ray, epsilon = 1e-10) -> np.ndarray[np.floating, Any]: + def vecdot(a:Vector3Array|Vector3, b:Vector3Array) -> FloatArray: return np.sum(a*b, axis=-1).view(FloatArray) # Translated from https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm#C++_implementation - triangle_array_comp = triangle_array.swapaxes(0, 1) - v0, v1, v2 = triangle_array_comp[0], triangle_array_comp[1], triangle_array_comp[2] - + v0, v1, v2 = triangle_array.verticies + # np.subtract(v1, v0, out= v1) # np.subtract(v2, v0, out= v2) v1 = v1 - v0 v2 = v2 - v0 - ray_cross_e2 = np.cross(ray.direction, v2) - det = np.vecdot(v1, ray_cross_e2) + ray_cross_e2 = np.cross(ray.direction, v2).view(Vector3Array) + det = vecdot(v1, ray_cross_e2) # if -epsilon < det < epsilon: # return None # This ray is parallel to this triangle. old_error_state = np.seterr(divide='ignore', invalid='ignore') inv_det = 1.0 / det - s = ray.origin - v0 - u = inv_det * np.vecdot(s, ray_cross_e2) + s = (ray.origin - v0).view(Vector3Array) + u = inv_det * vecdot(s, ray_cross_e2) - s_cross_e1 = np.cross(s, v1) - v = inv_det * np.vecdot(ray.direction, s_cross_e1) + s_cross_e1 = np.cross(s, v1).view(Vector3Array) + v = inv_det * vecdot(ray.direction, s_cross_e1) valid = u < -epsilon np.logical_or(valid, 1 < u, out=valid) @@ -108,14 +108,24 @@ def batch_triangle_ray_intersection(triangle_array, ray:Ray, epsilon = 1e-10) -> np.logical_or(valid, 1 < u + v, out=valid) np.invert(valid, out=valid) - t = np.empty(triangle_array.shape[0], dtype=np.float64) + t = np.empty(len(triangle_array), dtype=float) t.fill(-1) - np.multiply(inv_det, np.vecdot(v2, s_cross_e1), out = t, where = valid) + np.multiply(inv_det, vecdot(v2, s_cross_e1), out = t, where = valid) np.seterr(**old_error_state) return t + _centroids:Vector3Array|None = None + @property + def centroids(self) -> Vector3Array: + if self._centroids is None: + self._centroids = self.triangles.sum(axis = 1).view(Vector3Array) / 3 + return self._centroids + @centroids.setter + def centroids(self, value:Vector3Array) -> None: + self._centroids = value + def save(self, filename, override = False) -> None: with open(filename, 'wb' if override else 'xb') as file: pickle.dump(self, file) @@ -125,4 +135,4 @@ def load(filename) -> 'TriangleMesh': obj = pickle.load(file) if isinstance(obj, TriangleMesh): return obj - raise TypeError(f'The object saved in {filename} is type {type(obj)} not TriangleMesh') \ No newline at end of file + raise TypeError(f'The object saved in {filename} is type {type(obj)} not TriangleMesh') diff --git a/Raytrace/PlotRays.py b/Raytrace/PlotRays.py index 6fa0dd2..a28dccc 100644 --- a/Raytrace/PlotRays.py +++ b/Raytrace/PlotRays.py @@ -1,5 +1,5 @@ -from Raytrace.TriangleMesh import TriangleMesh -from Raytrace.SideScan import ScanReading, RawScanReading +from . import TriangleMesh +from .SideScan import ScanReading, RawScanReading import numpy as np try: import plotly.graph_objs as go # type: ignore diff --git a/Raytrace/SideScan.py b/Raytrace/SideScan.py index 83db3b6..1900053 100644 --- a/Raytrace/SideScan.py +++ b/Raytrace/SideScan.py @@ -1,4 +1,5 @@ -from Raytrace.TriangleMesh import TriangleMesh, Ray +from . import TriangleMesh +from CustomTypes import Ray, Vector3, Vector3Array from typing import Self import numpy as np import time @@ -34,7 +35,7 @@ def combine_min(self, other:Self) -> Self: new_distances[other_smaller] = other.distances[other_smaller] # HACK: new rays should be generated based on which distance was closer - new_rays = [Ray(i, j) for i, j in zip(self.directions, self.origins)] + new_rays = list(map(Ray, self.directions, self.origins)) return type(self)(new_distances, new_rays) def print_summary(self) -> None: @@ -85,7 +86,7 @@ def process_raw(self) -> None: ldist = norm - (np.arange(0,1,1/self.result_reselution) + 0.5/self.result_reselution).reshape((-1,1)) smooth_val = self.smooth_dist / (self.max_dist - self.min_dist) - lval = np.pow(np.maximum(0, np.square(smooth_val) - np.square(ldist)),3) / (32/35*smooth_val**7) / len(self.raw.distances) + lval = (np.maximum(0, np.square(smooth_val) - np.square(ldist))**3) / (32/35*smooth_val**7) / len(self.raw.distances) self.result = np.sum(lval, axis = 1) np.seterr(**old_error_state) @@ -130,13 +131,14 @@ def generate_rays(orientation:Ray, min_angle:float, max_angle:float, angle_resel angle_step = (max_angle-min_angle)/(angle_reselution-1) angles = np.arange(min_angle, max_angle + angle_step/2, angle_step) origin = orientation.origin - pitch = np.arctan2(orientation.direction[2],np.sqrt(np.dot(orientation.direction[0:2], orientation.direction[0:2]))) - yaw = np.arctan2(orientation.direction[1], orientation.direction[0]) # yaw = 0 => facing due +x + pitch:float = np.arctan2(orientation.direction[2],np.sqrt(np.dot(orientation.direction.nd[0:2], orientation.direction.nd[0:2]))) + yaw:float = np.arctan2(orientation.direction[1], orientation.direction[0]) # yaw = 0 => facing due +x # Precomputed rotation matrix sin_comp = np.array([-np.sin(pitch) * np.cos(yaw), np.sin(pitch) * np.sin(yaw), np.cos(pitch)]) cos_comp = np.array([np.sin(yaw), np.cos(yaw), 0]) - return [Ray(np.sin(a) * sin_comp + np.cos(a) * cos_comp, origin) for a in angles] + directions = Vector3Array(np.array([np.sin(a) * sin_comp + np.cos(a) * cos_comp for a in angles])) + return list(map(Ray, directions, [origin]*angle_reselution)) \ No newline at end of file diff --git a/Raytrace/__init__.py b/Raytrace/__init__.py index e69de29..9954358 100644 --- a/Raytrace/__init__.py +++ b/Raytrace/__init__.py @@ -0,0 +1,5 @@ +from .Meshs import ( + TriangleMesh as TriangleMesh, + BVHMesh as BVHMesh, + CompositeMesh as CompositeMesh, +) \ No newline at end of file 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..995f14f --- /dev/null +++ b/ShipGen/Generate.py @@ -0,0 +1,263 @@ +import os + +import numpy as np +from tqdm import tqdm +import torch +from stl import Mesh + +from CustomTypes import TriangleArray + +from .tools import Guided_Cond_DDPM_Tools as GC_DDPM +from .tools.HullParameterization import Hull_Parameterization as HP + +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) -> TriangleArray: + return mesh.vectors.view(TriangleArray) +def mesh_list_to_triangles(mesh_list:list[Mesh]) -> list[TriangleArray]: + 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]) + +