diff --git a/designs/OC4semi.yaml b/designs/OC4semi.yaml index c3068604..97533447 100644 --- a/designs/OC4semi.yaml +++ b/designs/OC4semi.yaml @@ -1060,10 +1060,13 @@ turbine: platform: + intersectMesh : 1 # [-] 0 for without and 1 for intersection min_freq_BEM : 0.03 # [Hz] lowest frequency and frequency interval to use in BEM analysis dz_BEM : 3.0 # [m] axial discretization panel length target for BEM analysis da_BEM : 2.0 # [m] azimuthal discretization panel length target for BEM analysis - potModMaster : 1 + potModMaster : 2 + characteristic_length_min: 0.5 # [m] Min mesh element size + characteristic_length_max: 0.9 # [m] Max mesh element size members: # list all members here @@ -1078,6 +1081,8 @@ platform: stations : [-20, 10 ] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : [ 6.5, 6.5] # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.03 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 0.0 # [m] length of extension on end A of the main column + extensionB : 0.0 # [m] length of extension on end B of the main column Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1102,6 +1107,8 @@ platform: stations : [-20, -14, -14, 12 ] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : [ 24, 24, 12, 12] # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.06 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 0.0 # [m] length of extension on both sides of the offset column + extensionB : 0.0 # [m] length of extension on both sides of the offset column Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1129,6 +1136,8 @@ platform: stations : [0, 1 ] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : [ 1.6, 1.6] # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.0175 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 5.0 # [m] length of extension on both sides of the pontoon + extensionB : 5.0 # [m] length of extension on both sides of the pontoon Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1148,6 +1157,8 @@ platform: stations : [0, 1 ] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : [ 1.6, 1.6] # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.0175 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 5.0 # [m] length of extension on both sides of the pontoon + extensionB : 5.0 # [m] length of extension on both sides of the pontoon Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1167,6 +1178,8 @@ platform: stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : [ 1.6, 1.6] # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.0175 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 5.0 # [m] length of extension on both sides of the pontoon + extensionB : 5.0 # [m] length of extension on both sides of the pontoon Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1186,6 +1199,8 @@ platform: stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : [ 1.6, 1.6] # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.0175 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 5.0 # [m] length of extension on both sides of the pontoon + extensionB : 5.0 # [m] length of extension on both sides of the pontoon Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1205,6 +1220,8 @@ platform: stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : [ 1.6, 1.6] # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.0175 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 3.0 # [m] length of extension on both sides of the pontoon + extensionB : 3.0 # [m] length of extension on both sides of the pontoon Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) diff --git a/designs/VolturnUS-S.yaml b/designs/VolturnUS-S.yaml index f7e543a0..0fc12331 100644 --- a/designs/VolturnUS-S.yaml +++ b/designs/VolturnUS-S.yaml @@ -1085,8 +1085,11 @@ turbine: platform: - potModMaster : 1 # [int] master switch for potMod variables; 0=keeps all member potMod vars the same, 1=turns all potMod vars to False (no HAMS), 2=turns all potMod vars to True (no strip) + intersectMesh : 0 # 0 for without and 1 for intersection + potModMaster : 0 # [int] master switch for potMod variables; 0=keeps all member potMod vars the same, 1=turns all potMod vars to False (no HAMS), 2=turns all potMod vars to True (no strip) dlsMax : 5.0 # maximum node splitting section amount for platform members; can't be 0 + characteristic_length_min: 0.5 # [m] Min mesh element size + characteristic_length_max: 0.9 # [m] Max mesh element size members: # list all members here @@ -1101,6 +1104,8 @@ platform: stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : 10.0 # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 0.0 # [m] length of extension on end A of the center column. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + extensionB : 0.0 # [m] length of extension on end B of the center column. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1124,6 +1129,8 @@ platform: stations : [0, 35] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : 12.5 # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 0.0 # [m] length of extension on end A of the outer column. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + extensionB : 0.0 # [m] length of extension on end B of the outer column. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1143,13 +1150,15 @@ platform: rA : [ 5 , 0, -16.5] # [m] end A coordinates rB : [ 45.5, 0, -16.5] # [m] and B coordinates heading : [ 60, 180, 300] # [deg] heading rotation of column about z axis (for repeated members) - shape : rect # [-] circular or rectangular + shape : rec # [-] circular or rectangular gamma : 0.0 # [deg] twist angle about the member's z-axis potMod : False # [bool] Whether to model the member with potential flow (BEM model) plus viscous drag or purely strip theory # --- outer shell including hydro--- stations : [0, 40.5] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : [12.5, 7.0] # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 5.0 # [m] length of extension on end A of the pontoon. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + extensionB : 5.0 # [m] length of extension on end B of the pontoon. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) @@ -1171,6 +1180,8 @@ platform: stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB d : 0.91 # [m] diameters if circular or side lengths if rectangular (can be pairs) t : 0.01 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 5.0 # [m] length of extension on end A of the pontoon + extensionB : 5.0 # [m] length of extension on end B of the pontoon Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) diff --git a/docs/usage.rst b/docs/usage.rst index 43bc6c89..abdfb830 100644 --- a/docs/usage.rst +++ b/docs/usage.rst @@ -286,9 +286,9 @@ Platform .. code-block:: python platform: - - potModMaster : 1 # [int] master switch for potMod variables; 0=keeps all member potMod vars the same, 1=turns all potMod vars to False (no HAMS), 2=turns all potMod vars to True (no strip) - dlsMax : 5.0 # maximum node splitting section amount for platform members; can't be 0 + intersectMesh : 1 # [int] 0 for disabling and 1 for enabling meshing for intersected members, make sure you install pygmsh and meshmagick before using this option + potModMaster : 1 # [int] master switch for potMod variables; 0=keeps all member potMod vars the same, 1=turns all potMod vars to False (no HAMS), 2=turns all potMod vars to True (no strip) + dlsMax : 5.0 # maximum node splitting section amount for platform members; can't be 0 members: # list all members here @@ -300,14 +300,16 @@ Platform gamma : 0.0 # [deg] twist angle about the member's z-axis potMod : True # [bool] Whether to model the member with potential flow (BEM model) plus viscous drag or purely strip theory # --- outer shell including hydro--- - stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB - d : 10.0 # [m] diameters if circular or side lengths if rectangular (can be pairs) - t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) - Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) - Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) - CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) - CaEnd : 0.6 # [-] end axial added mass coefficient (optional, scalar or list of same length as stations) - rho_shell : 7850 # [kg/m3] + stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB + d : 10.0 # [m] diameters if circular or side lengths if rectangular (can be pairs) + t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 0.0 # [m] length of extension on end A of the center column. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + extensionB : 0.0 # [m] length of extension on end B of the center column. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) + Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) + CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) + CaEnd : 0.6 # [-] end axial added mass coefficient (optional, scalar or list of same length as stations) + rho_shell : 7850 # [kg/m3] # --- handling of end caps or any internal structures if we need them --- cap_stations : [ 0 ] # [m] location along member of any inner structures (in same scaling as set by 'stations') cap_t : [ 0.001 ] # [m] thickness of any internal structures @@ -322,14 +324,16 @@ Platform gamma : 0.0 # [deg] twist angle about the member's z-axis potMod : True # [bool] Whether to model the member with potential flow (BEM model) plus viscous drag or purely strip theory # --- outer shell including hydro--- - stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB - d : 12.5 # [m] diameters if circular or side lengths if rectangular (can be pairs) - t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) - Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) - Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) - CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) - CaEnd : 0.6 # [-] end axial added mass coefficient (optional, scalar or list of same length as stations) - rho_shell : 7850 # [kg/m3] + stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB + d : 12.5 # [m] diameters if circular or side lengths if rectangular (can be pairs) + t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 0.0 # [m] length of extension on end A of the outer column. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + extensionB : 0.0 # [m] length of extension on end B of the outer column. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) + Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) + CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) + CaEnd : 0.6 # [-] end axial added mass coefficient (optional, scalar or list of same length as stations) + rho_shell : 7850 # [kg/m3] # --- ballast --- l_fill : 1.4 # [m] rho_fill : 5000 # [kg/m3] @@ -347,16 +351,18 @@ Platform gamma : 0.0 # [deg] twist angle about the member's z-axis potMod : False # [bool] Whether to model the member with potential flow (BEM model) plus viscous drag or purely strip theory # --- outer shell including hydro--- - stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB - d : [12.5, 7.0] # [m] diameters if circular or side lengths if rectangular (can be pairs) - t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) - Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) - Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) - CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) - CaEnd : 0.6 # [-] end axial added mass coefficient (optional, scalar or list of same length as stations) - rho_shell : 7850 # [kg/m3] - l_fill : 43.0 # [m] - rho_fill : 1025.0 # [kg/m3] + stations : [0, 1] # [-] location of stations along axis. Will be normalized such that start value maps to rA and end value to rB + d : [12.5, 7.0] # [m] diameters if circular or side lengths if rectangular (can be pairs) + t : 0.05 # [m] wall thicknesses (scalar or list of same length as stations) + extensionA : 5.0 # [m] length of extension on end A of the pontoon. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + extensionB : 5.0 # [m] length of extension on end B of the pontoon. This extension is to ensure a valid boolean union operation when members are intersected. This will be autonatically determined when using within WEIS. + Cd : 0.8 # [-] transverse drag coefficient (optional, scalar or list of same length as stations) + Ca : 1.0 # [-] transverse added mass coefficient (optional, scalar or list of same length as stations) + CdEnd : 0.6 # [-] end axial drag coefficient (optional, scalar or list of same length as stations) + CaEnd : 0.6 # [-] end axial added mass coefficient (optional, scalar or list of same length as stations) + rho_shell : 7850 # [kg/m3] + l_fill : 43.0 # [m] + rho_fill : 1025.0 # [kg/m3] - ... diff --git a/raft/IntersectionMesh.py b/raft/IntersectionMesh.py new file mode 100644 index 00000000..9ca16196 --- /dev/null +++ b/raft/IntersectionMesh.py @@ -0,0 +1,218 @@ +import math +import pygmsh +import meshmagick +import subprocess +import numpy as np +import os +import platform +import yaml + + +def load_yaml(file_path): + try: + with open(file_path, "r") as file: + yaml_data = yaml.safe_load(file) + return yaml_data + except Exception as e: + print(f"Error loading YAML file: {e}") + return None + + +def meshMember(geom, headings, rA, rB, radius, member_id=0, + stations=[0.0, 1.0], diameters=None, extensionA=0, extensionB=0): + cylinders = [] + + # Compute full axis and direction unit vector + axis_vec = [rB[i] - rA[i] for i in range(3)] + axis_length = math.sqrt(sum(a**2 for a in axis_vec)) + direction = [a / axis_length for a in axis_vec] + + # Apply extensions + rA_ext = [rA[i] - extensionA * direction[i] for i in range(3)] + rB_ext = [rB[i] + extensionB * direction[i] for i in range(3)] + axis_full = [rB_ext[i] - rA_ext[i] for i in range(3)] + + uniform = isinstance(diameters, (int, float)) or (isinstance(diameters, list) and len(diameters) == 1) + + for idx in range(len(headings)): + if np.all(diameters==diameters[0]): + start = rA_ext + end = rB_ext + axis_segment = [end[i] - start[i] for i in range(3)] + cone = geom.add_cylinder(start, axis_segment, diameters[0]/2) + label = f"Cylinder_{member_id}_{idx}" + print(f"Meshing {label} | Start: {start} | End: {end} | Radius: {diameters[0]/2}->{diameters[0]/2}") + + geom.add_physical(cone, label=label) + cylinders.append(cone) + else: + for s in range(len(stations) - 1): + t0 = (stations[s] - stations[0]) / (stations[-1] - stations[0]) + t1 = (stations[s + 1] - stations[0]) / (stations[-1] - stations[0]) + + if abs(t1 - t0) < 1e-6: + continue + + # ⏱ Interpolate segment along extended axis + start = [rA_ext[i] + t0 * axis_full[i] for i in range(3)] + end = [rA_ext[i] + t1 * axis_full[i] for i in range(3)] + axis_segment = [end[i] - start[i] for i in range(3)] + + if uniform or diameters is None: + radius_start = radius_end = radius + else: + radius_start = diameters[s] / 2 + radius_end = diameters[s + 1] / 2 + + label = f"Cylinder_{member_id}_{idx}_seg{s}" + print(f"Meshing {label} | Start: {start} | End: {end} | Radius: {radius_start}->{radius_end}") + + if abs(radius_start - radius_end) < 1e-6: + cone = geom.add_cylinder(start, axis_segment, radius_start) + else: + cone = geom.add_cone(start, axis_segment, radius_start, radius_end) + + geom.add_physical(cone, label=label) + cylinders.append(cone) + + return cylinders + + + +def meshRectangularMember(geom, heading, rA, rB, widths, heights, member_id=0, + stations=[0.0, 1.0], extensionA=0, extensionB=0): + boxes = [] + + # Compute axis and unit vector + axis_vec = [rB[i] - rA[i] for i in range(3)] + axis_length = math.sqrt(sum(a ** 2 for a in axis_vec)) + direction = [a / axis_length for a in axis_vec] + + # Apply extensions to both ends + rA_ext = [rA[i] - extensionA * direction[i] for i in range(3)] + rB_ext = [rB[i] + extensionB * direction[i] for i in range(3)] + axis_full = [rB_ext[i] - rA_ext[i] for i in range(3)] + + for s in range(len(stations) - 1): + t0 = (stations[s] - stations[0]) / (stations[-1] - stations[0]) + t1 = (stations[s + 1] - stations[0]) / (stations[-1] - stations[0]) + + if abs(t1 - t0) < 1e-6: + continue + + start = [rA_ext[i] + t0 * axis_full[i] for i in range(3)] + end = [rA_ext[i] + t1 * axis_full[i] for i in range(3)] + axis_segment = [end[i] - start[i] for i in range(3)] + + width = widths[s] + height = heights[s] + length = math.sqrt(sum(a ** 2 for a in axis_segment)) + box_size = [length, width, height] + + # Box at origin + box = geom.add_box( + [-box_size[0] / 2, -box_size[1] / 2, -box_size[2] / 2], + box_size + ) + + # Rotation to align + dx, dy = axis_segment[0], axis_segment[1] + theta = math.atan2(dy, dx) + geom.rotate(box, point=(0, 0, 0), angle=theta, axis=(0, 0, 1)) + + # Translate box to its center + box_center = [ + start[0] + axis_segment[0] / 2, + start[1] + axis_segment[1] / 2, + start[2] + axis_segment[2] / 2 + ] + geom.translate(box, box_center) + + label = f"Box_{member_id}_seg{s}" + print(f"Box: {label} | Center: {box_center} | Heading: {math.degrees(theta):.1f} | theta: {theta:.1f} | Size: {box_size}") + geom.add_physical(box, label=label) + boxes.append(box) + + return boxes + + +def mesh(meshDir=os.path.join(os.getcwd(),'BEM'), cylindrical_members=[], rectangular_members=[], dmin=0.1, dmax=1): + + print(f"Total cylindrical members: {len(cylindrical_members)}") + print(f"Total rectangular members: {len(rectangular_members)}") + all_shapes = [] + + with pygmsh.occ.Geometry() as geom: + geom.characteristic_length_min = dmin + geom.characteristic_length_max = dmax + + for member_id, cyl in enumerate(cylindrical_members): + try: + cylinders = meshMember( + geom, cyl.get("heading", [0]), cyl["rA"], cyl["rB"], cyl.get("radius"), + member_id=member_id, + stations=cyl.get("stations", [0.0, 1.0]), diameters=cyl.get("diameters"), extensionA=cyl.get("extensionA", 0), extensionB=cyl.get("extensionB", 0) + ) + all_shapes.extend(cylinders) + except Exception as e: + print(f"Failed to mesh cylindrical member {member_id}: {e}") + + for rect_id, rect in enumerate(rectangular_members): + try: + boxes = meshRectangularMember( + geom, rect.get("heading", [0]), rect["rA"], rect["rB"], + rect["widths"], rect["heights"], + member_id=rect_id, stations=rect.get("stations", [0.0, 1.0]), extensionA=rect.get("extensionA", 0), extensionB=rect.get("extensionB", 0) + ) + all_shapes.extend(boxes) + except Exception as e: + print(f"Failed to mesh rectangular member {rect_id}: {e}") + + if not all_shapes: + print("No shapes were added! Check your YAML or input logic.") + return + + if os.path.isdir(meshDir) is not True: + os.makedirs(meshDir) + + try: + combined = geom.boolean_union(all_shapes) + geom.add_physical(combined, label="CombinedGeometry") + mesh = geom.generate_mesh() + stl_path = os.path.join(meshDir, "Platform.stl") + mesh.write(stl_path) + except Exception as e: + print(f"Boolean union or meshing failed: {e}") + return + + try: + mesh_path = os.path.join(meshDir, "HullMesh.pnl") + intermediate_path = os.path.join(meshDir, "Platform.pnl") + if platform.system() == "Windows": + + subprocess.run([ + "meshmagick.exe", + stl_path, "-o", intermediate_path, + "--input-format", "stl", "--output-format", "pnl" + ], check=True) + + subprocess.run([ + "meshmagick.exe", + intermediate_path, "-c", "Oxy", "-o", mesh_path + ], check=True) + else: + subprocess.run([ + "meshmagick", + stl_path, "-o", intermediate_path, + "--input-format", "stl", "--output-format", "pnl" + ], check=True) + + subprocess.run([ + "meshmagick", + intermediate_path, "-c", "Oxy", "-o", mesh_path + ], check=True) + except subprocess.CalledProcessError as e: + print(f"Meshmagick failed: {e}") + +if __name__ == "__main__": + mesh() diff --git a/raft/member2pnl.py b/raft/member2pnl.py index 4f56f242..e1fd0304 100644 --- a/raft/member2pnl.py +++ b/raft/member2pnl.py @@ -499,6 +499,351 @@ def meshMemberForGDF(stations, diameters, rA, rB, dz_max=0, da_max=0, endA=True, return vertices2.T +# For rectangular members + +def meshRectangularMember(stations, widths, heights, rA, rB, dz_max=0, dw_max=0, dh_max=0, savedNodes=[], savedPanels=[]): + ''' + Creates mesh for a rectangular member. + + Parameters + ---------- + stations: list + locations along member axis at which the cross section will be specified + widths: list + corresponding widths along member + heights: list + corresponding heights along member + rA, rB: list + member end point coordinates + dz_max: float + maximum panel height + dw_max: float + maximum panel width + dh_max: float + maximum panel height + savedNodes : list of lists + all the node coordinates already saved + savedPanels : list + the information for all the panels already saved: panel_number number_of_vertices Vertex1_ID Vertex2_ID Vertex3_ID (Vertex4_ID) + + Returns + ------- + nodes : list + list of node coordinates + panels : list + the information for all the panels: panel_number number_of_vertices Vertex1_ID Vertex2_ID Vertex3_ID (Vertex4_ID) + ''' + + # discretization defaults + if dz_max == 0: + dz_max = stations[-1] / 20 + if dw_max == 0: + dw_max = np.max(widths) / 8 + if dh_max == 0: + dh_max = np.max(heights) / 8 + + # calculate member rotation matrix + rAB = rB - rA # displacement vector from end A to end B [m] + beta = np.arctan2(rAB[1], rAB[0]) # member incline heading from x axis + phi = np.arctan2(np.sqrt(rAB[0]**2 + rAB[1]**2), rAB[2]) # member incline angle from vertical + + s1 = np.sin(beta) + c1 = np.cos(beta) + s2 = np.sin(phi) + c2 = np.cos(phi) + s3 = np.sin(0.0) + c3 = np.cos(0.0) + + R = np.array([[c1 * c2 * c3 - s1 * s3, -c3 * s1 - c1 * c2 * s3, c1 * s2], + [c1 * s3 + c2 * c3 * s1, c1 * c3 - c2 * s1 * s3, s1 * s2], + [-c3 * s2, s2 * s3, c2]]) # Z1Y2Z3 from https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix + + nSavedPanelsOld = len(savedPanels) + + # Create intermediate stations based on dz_max + new_stations = [stations[0]] + new_widths = [widths[0]] + new_heights = [heights[0]] + for i in range(len(stations) - 1): + z1, z2 = stations[i], stations[i + 1] + w1, w2 = widths[i], widths[i + 1] + h1, h2 = heights[i], heights[i + 1] + + dz = z2 - z1 + num_divisions = max(int(np.ceil(dz / dz_max)), 1) + z_divisions = np.linspace(z1, z2, num_divisions + 1) + for j in range(1, num_divisions + 1): + new_stations.append(z_divisions[j]) + new_widths.append(w1 + (w2 - w1) * (z_divisions[j] - z1) / dz) + new_heights.append(h1 + (h2 - h1) * (z_divisions[j] - z1) / dz) + + # Generate panels for all divisions + for i in range(len(new_stations) - 1): + w1, w2 = new_widths[i], new_widths[i + 1] + h1, h2 = new_heights[i], new_heights[i + 1] + z1, z2 = new_stations[i], new_stations[i + 1] + + nw = max(int(np.ceil(w1 / dw_max)), int(np.ceil(w2 / dw_max))) + nh = max(int(np.ceil(h1 / dh_max)), int(np.ceil(h2 / dh_max))) + + x_divisions_w1 = np.linspace(-w1 / 2, w1 / 2, nw + 1) + x_divisions_w2 = np.linspace(-w2 / 2, w2 / 2, nw + 1) + y_divisions_h1 = np.linspace(-h1 / 2, h1 / 2, nh + 1) + y_divisions_h2 = np.linspace(-h2 / 2, h2 / 2, nh + 1) + + for iw in range(nw): + for ih in range(nh): + # Panel 1: Front face + X = [x_divisions_w1[iw], x_divisions_w1[iw + 1], x_divisions_w2[iw + 1], x_divisions_w2[iw]] + Y = [-h1 / 2, -h1 / 2, -h2 / 2, -h2 / 2] + Z = [z1, z1, z2, z2] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + makePanel(nodes[0], nodes[1], nodes[2], savedNodes, savedPanels) + + # Panel 2: Back face + X = [x_divisions_w1[iw], x_divisions_w1[iw + 1], x_divisions_w2[iw + 1], x_divisions_w2[iw]] + Y = [h1 / 2, h1 / 2, h2 / 2, h2 / 2] + Z = [z1, z1, z2, z2] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + makePanel(nodes[0], nodes[1], nodes[2], savedNodes, savedPanels) + + for iw in range(nw): + for ih in range(nh): + # Panel 3: Left face + X = [-w1 / 2, -w1 / 2, -w2 / 2, -w2 / 2] + Y = [y_divisions_h1[ih], y_divisions_h1[ih + 1], y_divisions_h2[ih + 1], y_divisions_h2[ih]] + Z = [z1, z1, z2, z2] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + makePanel(nodes[0], nodes[1], nodes[2], savedNodes, savedPanels) + + # Panel 4: Right face + X = [w1 / 2, w1 / 2, w2 / 2, w2 / 2] + Y = [y_divisions_h1[ih], y_divisions_h1[ih + 1], y_divisions_h2[ih + 1], y_divisions_h2[ih]] + Z = [z1, z1, z2, z2] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + makePanel(nodes[0], nodes[1], nodes[2], savedNodes, savedPanels) + + # Mesh the end faces + wA = new_widths[0] + hA = new_heights[0] + wB = new_widths[-1] + hB = new_heights[-1] + zA = new_stations[0] + zB = new_stations[-1] + + # End A face + nwA = max(int(np.ceil(wA / dw_max)), 1) + nhA = max(int(np.ceil(hA / dh_max)), 1) + x_divisions_wA = np.linspace(-wA / 2, wA / 2, nwA + 1) + y_divisions_hA = np.linspace(-hA / 2, hA / 2, nhA + 1) + + for iw in range(nwA): + for ih in range(nhA): + X = [x_divisions_wA[iw], x_divisions_wA[iw + 1], x_divisions_wA[iw + 1], x_divisions_wA[iw]] + Y = [y_divisions_hA[ih], y_divisions_hA[ih], y_divisions_hA[ih + 1], y_divisions_hA[ih + 1]] + Z = [zA, zA, zA, zA] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + makePanel(nodes[0], nodes[1], nodes[2], savedNodes, savedPanels) + + # End B face + nwB = max(int(np.ceil(wB / dw_max)), 1) + nhB = max(int(np.ceil(hB / dh_max)), 1) + x_divisions_wB = np.linspace(-wB / 2, wB / 2, nwB + 1) + y_divisions_hB = np.linspace(-hB / 2, hB / 2, nhB + 1) + + for iw in range(nwB): + for ih in range(nhB): + X = [x_divisions_wB[iw], x_divisions_wB[iw + 1], x_divisions_wB[iw + 1], x_divisions_wB[iw]] + Y = [y_divisions_hB[ih], y_divisions_hB[ih], y_divisions_hB[ih + 1], y_divisions_hB[ih + 1]] + Z = [zB, zB, zB, zB] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + makePanel(nodes[0], nodes[1], nodes[2], savedNodes, savedPanels) + + print(f'Of {len(savedPanels) - nSavedPanelsOld} generated panels, {len(savedPanels) - nSavedPanelsOld} were submerged and have been used in the mesh.') + + return savedNodes, savedPanels + +def meshRectangularMemberForGDF(stations, widths, heights, rA, rB, dz_max=0, dw_max=0, dh_max=0, endA=True, endB=True): + ''' + Creates mesh for a rectangular member. + + Parameters + ---------- + stations: list of locations along member axis at which the cross section will be specified + widths: list of corresponding widths along member + heights: list of corresponding heights along member + rA, rB: member end point coordinates + dz_max: maximum panel height + dw_max: maximum panel width + dh_max: maximum panel height + endA/endB: flag for whether to mesh each end + + Returns + ------- + vertices : array + An array containing the mesh point coordinates, size [3, 4*npanel] + ''' + + if len(stations) != len(widths) or len(stations) != len(heights): + raise ValueError("The lengths of stations, widths, and heights must be the same.") + + if dz_max == 0: + dz_max = stations[-1] / 20 + if dw_max == 0: + dw_max = max(widths) / 8 + if dh_max == 0: + dh_max = max(heights) / 8 + + rAB = np.array(rB) - np.array(rA) + beta = np.arctan2(rAB[1], rAB[0]) + phi = np.arctan2(np.sqrt(rAB[0]**2 + rAB[1]**2), rAB[2]) + + s1 = np.sin(beta) + c1 = np.cos(beta) + s2 = np.sin(phi) + c2 = np.cos(phi) + s3 = np.sin(0.0) + c3 = np.cos(0.0) + + R = np.array([[c1 * c2 * c3 - s1 * s3, -c3 * s1 - c1 * c2 * s3, c1 * s2], + [c1 * s3 + c2 * c3 * s1, c1 * c3 - c2 * s1 * s3, s1 * s2], + [-c3 * s2, s2 * s3, c2]]) + + nSavedPanelsOld = 0 + + new_stations = [stations[0]] + new_widths = [widths[0]] + new_heights = [heights[0]] + for i in range(len(stations) - 1): + z1, z2 = stations[i], stations[i + 1] + w1, w2 = widths[i], widths[i + 1] + h1, h2 = heights[i], heights[i + 1] + + dz = z2 - z1 + if dz == 0: + continue + + num_divisions = max(int(np.ceil(dz / dz_max)), 1) + z_divisions = np.linspace(z1, z2, num_divisions + 1) + for j in range(1, num_divisions + 1): + new_stations.append(z_divisions[j]) + new_widths.append(w1 + (w2 - w1) * (z_divisions[j] - z1) / dz) + new_heights.append(h1 + (h2 - h1) * (z_divisions[j] - z1) / dz) + + x, y, z = [], [], [] + + for i in range(len(new_stations) - 1): + w1, w2 = new_widths[i], new_widths[i + 1] + h1, h2 = new_heights[i], new_heights[i + 1] + z1, z2 = new_stations[i], new_stations[i + 1] + + if w1 == 0 or w2 == 0 or h1 == 0 or h2 == 0: + continue + + if dw_max == 0 or dh_max == 0: + raise ValueError("dw_max and dh_max must be non-zero.") + + nw = max(int(np.ceil(w1 / dw_max)), int(np.ceil(w2 / dw_max))) + nh = max(int(np.ceil(h1 / dh_max)), int(np.ceil(h2 / dh_max))) + + x_divisions_w1 = np.linspace(-w1 / 2, w1 / 2, nw + 1) + x_divisions_w2 = np.linspace(-w2 / 2, w2 / 2, nw + 1) + y_divisions_h1 = np.linspace(-h1 / 2, h1 / 2, nh + 1) + y_divisions_h2 = np.linspace(-h2 / 2, h2 / 2, nh + 1) + + for iw in range(nw): + for ih in range(nh): + X = [x_divisions_w1[iw], x_divisions_w1[iw + 1], x_divisions_w2[iw + 1], x_divisions_w2[iw]] + Y = [-h1 / 2, -h1 / 2, -h2 / 2, -h2 / 2] + Z = [z1, z1, z2, z2] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + x.extend(nodes[0]) + y.extend(nodes[1]) + z.extend(nodes[2]) + + X = [x_divisions_w1[iw], x_divisions_w1[iw + 1], x_divisions_w2[iw + 1], x_divisions_w2[iw]] + Y = [h1 / 2, h1 / 2, h2 / 2, h2 / 2] + Z = [z1, z1, z2, z2] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + x.extend(nodes[0]) + y.extend(nodes[1]) + z.extend(nodes[2]) + + for iw in range(nw): + for ih in range(nh): + X = [-w1 / 2, -w1 / 2, -w2 / 2, -w2 / 2] + Y = [y_divisions_h1[ih], y_divisions_h1[ih + 1], y_divisions_h2[ih + 1], y_divisions_h2[ih]] + Z = [z1, z1, z2, z2] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + x.extend(nodes[0]) + y.extend(nodes[1]) + z.extend(nodes[2]) + + X = [w1 / 2, w1 / 2, w2 / 2, w2 / 2] + Y = [y_divisions_h1[ih], y_divisions_h1[ih + 1], y_divisions_h2[ih + 1], y_divisions_h2[ih]] + Z = [z1, z1, z2, z2] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + x.extend(nodes[0]) + y.extend(nodes[1]) + z.extend(nodes[2]) + + wA = new_widths[0] + hA = new_heights[0] + wB = new_widths[-1] + hB = new_heights[-1] + zA = new_stations[0] + zB = new_stations[-1] + + if endA and wA > 0 and hA > 0: + nwA = max(int(np.ceil(wA / dw_max)), 1) + nhA = max(int(np.ceil(hA / dh_max)), 1) + x_divisions_wA = np.linspace(-wA / 2, wA / 2, nwA + 1) + y_divisions_hA = np.linspace(-hA / 2, hA / 2, nhA + 1) + + for iw in range(nwA): + for ih in range(nhA): + X = [x_divisions_wA[iw], x_divisions_wA[iw + 1], x_divisions_wA[iw + 1], x_divisions_wA[iw]] + Y = [y_divisions_hA[ih], y_divisions_hA[ih], y_divisions_hA[ih + 1], y_divisions_hA[ih + 1]] + Z = [zA, zA, zA, zA] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + x.extend(nodes[0]) + y.extend(nodes[1]) + z.extend(nodes[2]) + + if endB and wB > 0 and hB > 0: + nwB = max(int(np.ceil(wB / dw_max)), 1) + nhB = max(int(np.ceil(hB / dh_max)), 1) + x_divisions_wB = np.linspace(-wB / 2, wB / 2, nwB + 1) + y_divisions_hB = np.linspace(-hB / 2, hB / 2, nhB + 1) + + for iw in range(nwB): + for ih in range(nhB): + X = [x_divisions_wB[iw], x_divisions_wB[iw + 1], x_divisions_wB[iw + 1], x_divisions_wB[iw]] + Y = [y_divisions_hB[ih], y_divisions_hB[ih], y_divisions_hB[ih + 1], y_divisions_hB[ih + 1]] + Z = [zB, zB, zB, zB] + nodes0 = np.array([X, Y, Z]) + nodes = np.matmul(R, nodes0) + rA[:, None] + x.extend(nodes[0]) + y.extend(nodes[1]) + z.extend(nodes[2]) + + vertices = np.array([x, y, z]) + + vertices2 = np.matmul(R, vertices) + rA[:, None] + + return vertices2.T + def writeMeshToGDF(vertices, filename="platform.gdf", aboveWater=True): npan = int(vertices.shape[0]/4) @@ -532,15 +877,33 @@ def writeMeshToGDF(vertices, filename="platform.gdf", aboveWater=True): if __name__ == "__main__": + + shape = "rectangular" + + if shape == "circular": - stations = [0, 6, 6,32] - diameters = [24, 24, 12, 12] - - rA = np.array([0, 0,-20]) - rB = np.array([0, 0, 12]) - - #nodes, panels = meshMember(stations, diameters, rA, rB, dz_max=1, da_max=1) - #writeMesh(nodes, panels) - - vertices = meshMemberForGDF(stations, diameters, rA, rB, dz_max=5, da_max=5) - writeMeshToGDF(vertices) + stations = [0, 6, 6,32] + diameters = [24, 24, 12, 12] + + rA = np.array([0, 0,-20]) + rB = np.array([0, 0, 12]) + + #nodes, panels = meshMember(stations, diameters, rA, rB, dz_max=1, da_max=1) + #writeMesh(nodes, panels) + + vertices = meshMemberForGDF(stations, diameters, rA, rB, dz_max=5, da_max=5) + writeMeshToGDF(vertices) + + elif shape == "rectangular": + stations= [-120, -12, -4, 10] + widths = [24, 24, 24, 24] + heights = [24, 24, 24, 24] + + rA = np.array([0, 0,-10]) + rB = np.array([0, 0, 22]) + + + vertices = meshRectangularMemberForGDF(stations, widths, heights, rA, rB, dz_max=20, dw_max=20, dh_max=20) + writeMeshToGDF(vertices) + savedNodes, savedPanels = meshRectangularMember(stations, widths, heights, rA, rB, dz_max=20, dw_max=20, dh_max=20) + writeMesh(savedNodes, savedPanels, "Test") diff --git a/raft/omdao_raft.py b/raft/omdao_raft.py index 3345faa5..a488b125 100644 --- a/raft/omdao_raft.py +++ b/raft/omdao_raft.py @@ -354,6 +354,8 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): member_scalar_t = members_opt['scalar_thicknesses'] member_scalar_d = members_opt['scalar_diameters'] member_scalar_coeff = members_opt['scalar_coefficients'] + intersectMesh = modeling_opt['intersection_mesh'] + nlines = mooring_opt['nlines'] nline_types = mooring_opt['nline_types'] @@ -501,6 +503,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): design['platform'] = {} design['platform']['potModMaster'] = int(modeling_opt['potential_model_override']) design['platform']['dlsMax'] = float(modeling_opt['dls_max']) + design['platform']["intersectMesh"] = intersectMesh # lowest BEM freq needs to be just below RAFT min_freq because of interpolation in RAFT if float(modeling_opt['min_freq_BEM']) >= modeling_opt['min_freq']: modeling_opt['min_freq_BEM'] = modeling_opt['min_freq'] - 1e-7 @@ -525,10 +528,14 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): mnpts = len(idx) rA = rA_0 + s_ghostA*(rB_0-rA_0) rB = rA_0 + s_ghostB*(rB_0-rA_0) + extensionA = np.linalg.norm(rA_0-rA) + extensionB = np.linalg.norm(rB_0-rB) design['platform']['members'][i]['name'] = m_name design['platform']['members'][i]['type'] = i + 2 design['platform']['members'][i]['rA'] = rA design['platform']['members'][i]['rB'] = rB + design['platform']['members'][i]['extensionA'] = extensionA + design['platform']['members'][i]['extensionB'] = extensionB design['platform']['members'][i]['shape'] = m_shape design['platform']['members'][i]['gamma'] = float(inputs[m_name+'gamma'][0]) design['platform']['members'][i]['potMod'] = members_opt[m_name+'potMod'] @@ -743,7 +750,33 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): # option to run level 1 load cases if True: #processCases: model.analyzeCases(meshDir=modeling_opt['BEM_dir']) + + # Plot mesh + if modeling_opt['plot_designs'] and intersectMesh: + try: + import trimesh + import matplotlib.pyplot as plt + except: + raise ImportError('trimesh not installed, please install trimesh to plot mesh') + mesh = trimesh.load(os.path.join(modeling_opt['BEM_dir'],"Input/Platform.stl")) + # Create a figure and axes + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + + # Get mesh vertices and faces + verts = mesh.vertices + faces = mesh.faces + + # Plot the mesh + ax.plot_trisurf( + verts[:, 0], verts[:, 1], faces, verts[:, 2], + color='lightblue', edgecolor='black', linewidth=0.2, alpha=0.8 + ) + ax.set_box_aspect([ub - lb for lb, ub in (getattr(ax, f'get_{a}lim')() for a in 'xyz')]) + plt.savefig(os.path.join(modeling_opt['BEM_dir'],"mesh_plot.png"), pad_inches=0, dpi=300) + + # get and process results results = model.calcOutputs() # Pattern matching for "responses" annd "properties" diff --git a/raft/raft_fowt.py b/raft/raft_fowt.py index d7a59bcf..a8ff9c09 100644 --- a/raft/raft_fowt.py +++ b/raft/raft_fowt.py @@ -12,6 +12,25 @@ import moorpy as mp from moorpy.helpers import lines2ss +# Attempt to import pygmsh and meshmagick with warnings if not installed +try: + import pygmsh +except ImportError: + pygmsh = None + print("Warning: 'pygmsh' is not installed. Meshing intersected members will not be available. Install it using 'pip install pygmsh==7.1.17'.") + +try: + import meshmagick +except ImportError: + meshmagick = None + print("Warning: 'meshmagick' is not installed. Meshing intersected members will not be available. Install it using 'pip install https://github.com/LHEEA/meshmagick/archive/refs/tags/3.4.zip'.") + +try: + import trimesh +except ImportError: + trimesh = None + print("Warning: 'trimesh' is not installed. Automatically plotting intermediate BEM meshes will not be available. Install it using 'pip install trimesh'.") + # deleted call to ccblade in this file, since it is called in raft_rotor # also ignoring changes to solveEquilibrium3 in raft_model and the re-addition of n=len(stations) in raft_member, based on raft_patch @@ -52,6 +71,10 @@ def __init__(self, design, w, mpb, depth=600, x_ref=0, y_ref=0, heading_adjust=0 self.Xi0 = np.zeros( self.nDOF) # mean offsets of platform from its reference point [m, rad] self.Xi = np.zeros([self.nDOF, self.nw], dtype=complex) # complex response amplitudes as a function of frequency [m, rad] self.heading_adjust = heading_adjust # rotation to the heading of the platform and mooring system to be applied [deg] + self.design = design + self.characteristic_length_min = design['platform'].get('characteristic_length_min', 1) + self.characteristic_length_max = design['platform'].get('characteristic_length_max', 3) + print(f"Characteristic lengths: min={self.characteristic_length_min}, max={self.characteristic_length_max}") # position in the array self.x_ref = x_ref # reference x position of the FOWT in the array [m] @@ -574,7 +597,7 @@ def calcStatics(self): self.props['Izz_sub'] = M_sub[5,5] - def calcBEM(self, dw=0, wMax=0, wInf=10.0, dz=0, da=0, headings=[0], meshDir=os.path.join(os.getcwd(),'BEM')): + def calcBEM(self, dw=0, wMax=0, wInf=10.0, dz=0, da=0, dh=0, headings=[0], meshDir=os.path.join(os.getcwd(),'BEM'), dmin=1, dmax=3): '''This generates a mesh for the platform and runs a BEM analysis on it using pyHAMS. It can also write adjusted .1 and .3 output files suitable for use with OpenFAST. @@ -612,23 +635,85 @@ def calcBEM(self, dw=0, wMax=0, wInf=10.0, dz=0, da=0, headings=[0], meshDir=os. dz = self.dz_BEM if dz==0 else dz # allow override if provided da = self.da_BEM if da==0 else da + #dh = self.dh_BEM if dh==0 else dh + if self.design["platform"]["intersectMesh"] == 0: + for mem in self.memberList: + if mem.potMod: # >>>>>>>>>>>>>>>> now using for rectnagular member and the dimensions are hardcoded. need to integrate with the .yaml file input + if mem.shape == "circular": + pnl.meshMember(mem.stations, mem.d, mem.rA, mem.rB, dz_max=dz, da_max=da, savedNodes=nodes, savedPanels=panels) + # for GDF output + vertices_i = pnl.meshMemberForGDF(mem.stations, mem.d, mem.rA, mem.rB, dz_max=dz, da_max=da) + vertices = np.vstack([vertices, vertices_i]) # append the member's vertices to the master list + elif mem.shape == "rectangular": + widths = mem.sl[:, 0] + heights = mem.sl[:, 1] + pnl.meshRectangularMember(mem.stations, widths, heights, mem.rA, mem.rB, dz_max=dz, dw_max=da, dh_max=dh, savedNodes=nodes, savedPanels=panels) # + # for GDF output + vertices_i = pnl.meshRectangularMemberForGDF(mem.stations, widths, heights, mem.rA, mem.rB, dz_max=dz, dw_max=da, dh_max=dh) + vertices = np.vstack([vertices, vertices_i]) # append the member's vertices to the master list + - for mem in self.memberList: - if mem.potMod: - pnl.meshMember(mem.stations, mem.d, mem.rA, mem.rB, - dz_max=dz, da_max=da, savedNodes=nodes, savedPanels=panels) + elif self.design["platform"]["intersectMesh"] == 1: + if pygmsh is None or meshmagick is None: + raise ImportError("The 'intersectMesh' option requires 'pygmsh' and 'meshmagick'. Please install them.") + + import raft.IntersectionMesh as intersectMesh + + cylindrical_members = [] + rectangular_members = [] - # for GDF output - vertices_i = pnl.meshMemberForGDF(mem.stations, mem.d, mem.rA, mem.rB, dz_max=dz, da_max=da) - vertices = np.vstack([vertices, vertices_i]) # append the member's vertices to the master list + platform = self.design.get("platform", {}) + for i, mem in enumerate(self.memberList): + if not mem.potMod: + continue + + stations = mem.stations + rA = mem.rA_original if hasattr(mem, "rA_original") else mem.rA + rB = mem.rB_original if hasattr(mem, "rB_original") else mem.rB + headings = mem.heading if hasattr(mem, "heading") else [0] + #print("name from raft:", mem.name) + #print("rA from raft: ",mem.rA) + #print("rB from raft: ", mem.rB) + #print(headings) + if mem.shape == "circular": + radius = mem.d / 2 if isinstance(mem.d, (int, float)) else None + diameters = mem.d + extensionA = getattr(mem, "extensionA", 0) + extensionB = getattr(mem, "extensionB", 0) + cylindrical_members.append({ + "rA": rA, + "rB": rB, + "radius": radius, + "heading": headings, + "stations": stations, + "diameters": diameters, + "extensionA": extensionA, + "extensionB": extensionB + }) + + elif mem.shape == "rectangular": + widths = mem.sl[:, 0].tolist() + heights = mem.sl[:, 1].tolist() + extensionA = getattr(mem, "extensionA", 0) + extensionB = getattr(mem, "extensionB", 0) + #print(f"Extension for {mem.name}: {extensionA}") + rectangular_members.append({ + "rA": rA, + "rB": rB, + "widths": widths, + "heights": heights, + "heading": headings, + "stations": stations, + "extensionA": extensionA, + "extensionB": extensionB + }) + + + intersectMesh.mesh(meshDir=os.path.join(meshDir,'Input'), cylindrical_members=cylindrical_members, rectangular_members=rectangular_members, dmin=self.characteristic_length_min, dmax=self.characteristic_length_max) if len(panels) == 0: print("WARNING: no panels to mesh.") - pnl.writeMesh(nodes, panels, oDir=os.path.join(meshDir,'Input')) # generate a mesh file in the HAMS .pnl format - - #pnl.writeMeshToGDF(vertices) # also a GDF for visualization - ph.create_hams_dirs(meshDir) # # HAMS needs a hydrostatics file, it's unused for .1 and .3, @@ -666,7 +751,6 @@ def calcBEM(self, dw=0, wMax=0, wInf=10.0, dz=0, da=0, headings=[0], meshDir=os. else: return - import pyhams.pyhams as ph # read the HAMS WAMIT-style output files addedMass, damping, w1 = ph.read_wamit1(hydroPath+'.1', TFlag=True) # first two entries in frequency dimension are expected to be zero-frequency then infinite frequency diff --git a/raft/raft_member.py b/raft/raft_member.py index 35c338c8..e9dcf94b 100644 --- a/raft/raft_member.py +++ b/raft/raft_member.py @@ -47,7 +47,8 @@ def __init__(self, mi, nw, BEM=[], heading=0): self.potMod = getFromDict(mi, 'potMod', dtype=bool, default=False) # hard coding BEM analysis enabled for now <<<< need to move this to the member YAML input instead <<< self.MCF = getFromDict(mi, 'MCF', dtype=bool, default=False) # Flag to use MacCamy-Fuchs correction or not - + self.extensionA = getFromDict(mi, 'extensionA', default=0) + self.extensionB = getFromDict(mi, 'extensionB', default=0) self.gamma = getFromDict(mi, 'gamma', default=0.) # twist angle about the member's z-axis [degrees] (if gamma=90, then the side lengths are flipped) rAB = self.rB0-self.rA0 # The relative coordinates of upper node from lower node [m] self.l = np.linalg.norm(rAB) # member length [m] diff --git a/raft/raft_model.py b/raft/raft_model.py index 4b776abb..4f4c97d3 100644 --- a/raft/raft_model.py +++ b/raft/raft_model.py @@ -17,7 +17,7 @@ import moorpy as mp import raft.raft_fowt as fowt from raft.helpers import * -from moorpy.helpers import dsolve2, set_axes_equal, dsolvePlot, lines2ss +from moorpy.helpers import dsolve2, set_axes_equal, dsolvePlot import copy #import F6T1RNA as structural # import turbine structural model functions diff --git a/tests/test_data/weis_options.yaml b/tests/test_data/weis_options.yaml index 5da9c03e..7888a538 100644 --- a/tests/test_data/weis_options.yaml +++ b/tests/test_data/weis_options.yaml @@ -15,6 +15,7 @@ modeling_options: runPyHAMS: true BEM_dir: /Users/dzalkind/Tools/WEIS-Main/examples/15_RAFT_Studies/none/BEM model_potential: [false, false, false, false, false, false, false, false, false, false] + intersection_mesh: 0 nfreq: 20 n_cases: 98 raft_dlcs: diff --git a/tests/test_omdao_OC3spar.py b/tests/test_omdao_OC3spar.py index b045295e..1298b08d 100644 --- a/tests/test_omdao_OC3spar.py +++ b/tests/test_omdao_OC3spar.py @@ -31,6 +31,7 @@ def testSpiderRegression(self): opt['modeling']['model_potential'] = [True,False, False,False] opt['modeling']['nfreq'] = 10 opt['modeling']['n_cases'] = 2 + opt['modeling']['intersection_mesh'] = 0 opt['turbine'] = {} opt['turbine']['npts'] = 11