13380864473

Committed 18 Feb 2025 12:38AM UTC coverage: 54.799% (-2.2%) from 57.037%

Build # 13380864473

Build Type

push

github

Committed by

web-flow

Commit Message

WEIS Input Visualization (#334)

* fix bug of args parse while running app

* delete unnecessary prints

* change horizontal subplots to vertical ones

* initial documentation

* delete readme file

* minor changes on graph layout

* update on kestrel set up

* merge weis viz docs into existing weis docs

* delete initial weis viz documentation

* update on docs after changing opt type setting

* build infrastructure for loading/viewing stl file

* add viz utils for windio - random mesh rendering

* initial attempt to load geometry variables and render cylinder

* update environment for geometry representation

* set interface to playaround meshes

* backup from hpc

* basic function to visualize single airfoil shape and polars

* allow multiple geom files and set the airfoil name with filename as dict keys

* add latex capabilities

* match color across plots and add symbols

* add re

* set same xy-ratio for airfoils shape and split airfoil polar plots

* add import page to upload weis input files

* update home page with card layout

* link airfoils loading func with import table

* add data preprocess function for 3d and make it accessible across pages

* change nickname into label

* init dev of windio blades

* import table update with deletion and handling error

* add LE/TE

* add tower page

* test tower cylinder with real geom parameters

* layout UI + render multiple towers from scratch

* move import page to home page

* update OpenFAST UI

* solve nonexistent object errors + dynamic card creation for OpenFast

* working blades

* debug nonexistent objects error + test with RAFT, OF opts

* polish layout and mathjax

* moving meshing to different file

* reconfig 3d layout

* show multiple comps multiple geoms combination at once

* show single global grids, removing duplicated cubedAxes

* add wt_options data to access whole data per file

* working 3d viz

* initial commit for viz github workflow

* test callbacks for 3d

* update workfl... (continued)

Run Details

0 of 493 new or added lines in 2 files covered. (0.0%)

4 existing lines in 1 file now uncovered.

6594 of 12033 relevant lines covered (54.8%)

0.55 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/weis/visualization/meshRender.py

NEW 1	import pandas as pd	×
NEW 2	import numpy as np	×
NEW 3	import matplotlib.pyplot as plt	×
NEW 4	from dash_vtk.utils import to_mesh_state	×
NEW 5	import pyvista as pv	×
6
NEW 7	try:	×
NEW 8	from vtkmodules.util.numpy_support import numpy_to_vtk, vtk_to_numpy	×
NEW 9	except:	×
NEW 10	from vtk.util.numpy_support import numpy_to_vtk, vtk_to_numpy	×
NEW 11	import yaml	×
12
13
NEW 14	def render_blade(turbineData, local=True):	×
15
16	# fetch airfoil names from array of airfoils
NEW 17	airfoils_by_names = {}	×
NEW 18	for a in turbineData['airfoils']:	×
NEW 19	airfoils_by_names[a['name']] = a	×
20
NEW 21	points = bladeMesh(turbineData['components']['blade'], airfoils_by_names)	×
22
NEW 23	if local:	×
NEW 24	mesh_state, mesh = render_our_own_delaunay(points)	×
25
26	else:
27	## transforming the blade to be in the global coordinate system
28
29	# Rotate about z-axis by 90-deg so that x is from HP to LP surface
NEW 30	points = rotation_transformation(points, [0, 0, np.pi/2])	×
31
32	# Translate to the hub position
NEW 33	points = translation_transformation(points, [0, 0, turbineData['components']['hub']['diameter'] / 2])	×
34
35	# if turbineData['assembly']['rotor_orientation'] has downwind (case insensitive) then rotate set downwindScalar to -1
NEW 36	downwindScalar = -1 if 'downwind' in turbineData['assembly']['rotor_orientation'].lower() else 1	×
37
38	# rotate about y-axis by the cone angle
NEW 39	coneAngle = turbineData['components']['hub']['cone_angle'] * -1 * downwindScalar # in rads is about x-axis	×
NEW 40	blade_1 = rotation_transformation(points, [0, coneAngle, 0])	×
41	# blade_2 = rotation_transformation(blade_1, [2*np.pi/3, 0, 0])
42	# blade_3 = rotation_transformation(blade_1, [4*np.pi/3, 0, 0])
43	# This achieves the symmetry of the rotor.
44
NEW 45	_, mesh_1 = render_our_own_delaunay(blade_1)	×
NEW 46	mesh_2 = mesh_1.rotate_x( 2np.pi/3 180/np.pi, point=(0,0,0), inplace=False)	×
NEW 47	mesh_3 = mesh_1.rotate_x( 4np.pi/3 180/np.pi, point=(0,0,0), inplace=False)	×
48
49
50	# angle by the tilt of the rotor
NEW 51	tiltAngle = turbineData['components']['nacelle']['drivetrain']['uptilt'] * downwindScalar # in rads is about y-axis	×
52	# blade_1 = rotation_transformation(blade_1, [0, tiltAngle, 0])
53	# blade_2 = rotation_transformation(blade_2, [0, tiltAngle, 0])
54	# blade_3 = rotation_transformation(blade_3, [0, tiltAngle, 0])
55
NEW 56	mesh_1 = mesh_1.rotate_y(tiltAngle * 180/np.pi, point=(0,0,0), inplace=False)	×
NEW 57	mesh_2 = mesh_2.rotate_y(tiltAngle * 180/np.pi, point=(0,0,0), inplace=False)	×
NEW 58	mesh_3 = mesh_3.rotate_y(tiltAngle * 180/np.pi, point=(0,0,0), inplace=False)	×
59
60	# Tanslation along the z-axis to the hub height
NEW 61	zTranslation = turbineData['assembly']['hub_height']	×
62	# Translation in the x-axis
NEW 63	xTranslation = turbineData['components']['nacelle']['drivetrain']['distance_tt_hub'] * downwindScalar * np.cos(tiltAngle) * -1	×
64
65	# blade_1 = translation_transformation(blade_1, [xTranslation, 0, zTranslation])
66	# blade_2 = translation_transformation(blade_2, [xTranslation, 0, zTranslation])
67	# blade_3 = translation_transformation(blade_3, [xTranslation, 0, zTranslation])
68
NEW 69	mesh_1 = mesh_1.translate((xTranslation, 0, zTranslation), inplace=False)	×
NEW 70	mesh_2 = mesh_2.translate((xTranslation, 0, zTranslation), inplace=False)	×
NEW 71	mesh_3 = mesh_3.translate((xTranslation, 0, zTranslation), inplace=False)	×
72
73	# _, mesh_1 = render_our_own_delaunay(blade_1)
74	# _, mesh_2 = render_our_own_delaunay(blade_2)
75	# _, mesh_3 = render_our_own_delaunay(blade_3)
76
NEW 77	mesh = mesh_1.merge(mesh_2).merge(mesh_3)	×
78
NEW 79	mesh_state = to_mesh_state(mesh)	×
80
81	# extracting all the points from the
NEW 82	rotor = mesh.points	×
83
NEW 84	extremes = extractExtremes(rotor)	×
85
NEW 86	return mesh_state, mesh, extremes	×
87
NEW 88	def render_hub(turbineData, local=True):	×
89
NEW 90	points = hubMesh(turbineData['components']['hub'])	×
91
NEW 92	if local:	×
NEW 93	mesh_state, mesh = render_our_own_delaunay(points)	×
94
95	else:
96	# if turbineData['assembly']['rotor_orientation'] has downwind (case insensitive) then rotate set downwindScalar to -1
NEW 97	downwindScalar = -1 if 'downwind' in turbineData['assembly']['rotor_orientation'].lower() else 1	×
98
99	# First rotate hub to align with x-axis (90 degrees about y),
100	# with additional 180 deg rotation for downwind configuration
NEW 101	points = rotation_transformation(points, [0, np.pi/2 + np.pi*int(downwindScalar==-1), 0])	×
102
103	# Apply tilt angle rotation
NEW 104	tiltAngle = turbineData['components']['nacelle']['drivetrain']['uptilt'] #* downwindScalar # in rads	×
NEW 105	points = rotation_transformation(points, [0, tiltAngle, 0])	×
106
107	# For downwind configuration, rotate 180 degrees about z-axis
NEW 108	if downwindScalar == -1:	×
NEW 109	points = rotation_transformation(points, [0, 0, np.pi])	×
110
111	# Tanslation along the z-axis to the hub height
NEW 112	zTranslation = turbineData['assembly']['hub_height']	×
113	# Translation in the x-axis (account for downwind configuration)
NEW 114	xTranslation = turbineData['components']['nacelle']['drivetrain']['distance_tt_hub'] * downwindScalar * np.cos(tiltAngle) * -1	×
115
NEW 116	points = translation_transformation(points, [xTranslation, 0, zTranslation])	×
117
NEW 118	mesh_state, mesh = render_our_own_delaunay(points)	×
119
NEW 120	extremes = extractExtremes(points)	×
121
NEW 122	return mesh_state, mesh, extremes	×
123
NEW 124	def render_Tower(turbineData, local=True):	×
125
NEW 126	points = towerMesh(turbineData['components']['tower'])	×
127
NEW 128	if local:	×
NEW 129	mesh_state, mesh = render_our_own_delaunay(points)	×
130	else:
131	# Tower does not have any orientation or tilt or translations
NEW 132	mesh_state, mesh = render_our_own_delaunay(points)	×
133
NEW 134	extremes = extractExtremes(points)	×
135
NEW 136	return mesh_state, mesh, extremes	×
137
138
NEW 139	def render_nacelle(turbineData, local=True):	×
140
NEW 141	points = nacelleMesh(turbineData['components']['nacelle'],turbineData['components']['hub'])	×
142
NEW 143	if local:	×
NEW 144	mesh_state, mesh = render_our_own_delaunay(points)	×
145	else:
146
147	# if turbineData['assembly']['rotor_orientation'] has downwind (case insensitive) then rotate set downwindScalar to -1
NEW 148	downwindScalar = -1 if 'downwind' in turbineData['assembly']['rotor_orientation'].lower() else 1	×
149
150	# need to translate it above the tower, the nacelle is always above the tower
NEW 151	zTranslation = turbineData['assembly']['hub_height']	×
152
153	# move the nacelle so that the -x face is touching the hub
NEW 154	xTranslation = turbineData['components']['nacelle']['drivetrain']['overhang'] * 0.5	×
155
156	# If downwind, rotate nacelle 180 degrees about z-axis
NEW 157	if downwindScalar == -1:	×
NEW 158	points = rotation_transformation(points, [0, 0, np.pi])	×
159
NEW 160	xTranslation = -xTranslation # Reverse the x-translation for downwind configuration	×
161
NEW 162	points = translation_transformation(np.array(points), [xTranslation, 0, zTranslation])	×
163
NEW 164	mesh_state, mesh = render_our_own_delaunay(points)	×
165
NEW 166	extremes = extractExtremes(points)	×
167
NEW 168	return mesh_state, mesh, extremes	×
169
NEW 170	def render_monopile(turbineData, local=True):	×
171
NEW 172	points = monopileMesh(turbineData['components']['monopile'])	×
173
NEW 174	if local:	×
NEW 175	mesh_state, mesh = render_our_own_delaunay(points)	×
176	else:
177	# Monopile does not have any orientation or tilt or translations
NEW 178	mesh_state, mesh = render_our_own_delaunay(points)	×
179
NEW 180	extremes = extractExtremes(points)	×
181
NEW 182	return mesh_state, mesh, extremes	×
183
NEW 184	def render_floatingPlatform(turbineData, local=True):	×
185
NEW 186	mesh = floatingPlatformMesh(turbineData['components']['floating_platform'],	×
187	turbineData['components']['mooring']) # this function unlike other components returns a mesh object
188
NEW 189	if local:	×
NEW 190	mesh_state = to_mesh_state(mesh)	×
191	else:
192	# Floating platform does not have any orientation or tilt or translations
NEW 193	mesh_state = to_mesh_state(mesh)	×
194
NEW 195	extremes = extractExtremes(mesh.points)	×
196
NEW 197	return mesh_state, mesh, extremes	×
198
199
NEW 200	def render_turbine(turbineData, components):	×
201
NEW 202	mesh_state_all = []	×
NEW 203	extremes_all = []	×
204
NEW 205	for idx, component in enumerate(components):	×
NEW 206	if component == 'blade':	×
NEW 207	mesh_state, mesh, extreme = render_blade(turbineData, local=False)	×
NEW 208	elif component == 'hub':	×
NEW 209	mesh_state, mesh, extreme = render_hub(turbineData, local=False)	×
NEW 210	elif component == 'tower':	×
NEW 211	mesh_state, mesh, extreme = render_Tower(turbineData, local=False)	×
NEW 212	elif component == 'nacelle':	×
NEW 213	mesh_state, mesh, extreme = render_nacelle(turbineData, local=False)	×
NEW 214	elif component == 'monopile':	×
NEW 215	mesh_state, mesh, extreme = render_monopile(turbineData, local=False)	×
NEW 216	elif component == 'floating_platform':	×
NEW 217	mesh_state, mesh, extreme = render_floatingPlatform(turbineData, local=False)	×
218	else:
NEW 219	raise ValueError(f"Component {component} not recognized")	×
220
NEW 221	mesh_state_all.append(mesh_state)	×
NEW 222	extremes_all.append(extreme)	×
NEW 223	if idx == 0:	×
NEW 224	mesh_all = mesh	×
225	else:
NEW 226	mesh_all = mesh_all.merge(mesh)	×
227
NEW 228	return mesh_state_all, mesh_all, extremes_all	×
229
230
231
NEW 232	def rotation_transformation(coords, angles):	×
233	'''
234	Apply rotation transformation to a set of coordinates
235	here, angle(1) is the rotation about x-axis, angle(2) is the rotation about y-axis, angle(3) is the rotation about z-axis
236	angles are in radians
237	'''
238	# Create rotation matrices for each angle
NEW 239	rot_x = np.array([[1, 0, 0], [0, np.cos(angles[0]), -np.sin(angles[0])], [0, np.sin(angles[0]), np.cos(angles[0])]])	×
NEW 240	rot_y = np.array([[np.cos(angles[1]), 0, np.sin(angles[1])], [0, 1, 0], [-np.sin(angles[1]), 0, np.cos(angles[1])]])	×
NEW 241	rot_z = np.array([[np.cos(angles[2]), -np.sin(angles[2]), 0], [np.sin(angles[2]), np.cos(angles[2]), 0], [0, 0, 1]])	×
242
243	# Apply rotation
NEW 244	coords = np.dot(rot_x, coords)	×
NEW 245	coords = np.dot(rot_y, coords)	×
NEW 246	coords = np.dot(rot_z, coords)	×
247
NEW 248	return coords	×
249
250
NEW 251	def translation_transformation(coords, translation):	×
252	'''
253	Apply translation transformation to a set of coordinates
254	'''
NEW 255	coords[0, :] += translation[0]	×
NEW 256	coords[1, :] += translation[1]	×
NEW 257	coords[2, :] += translation[2]	×
258
NEW 259	return coords	×
260
261
NEW 262	def extractExtremes(points):	×
263	'''
264	Extract the minimum and maximum values of the points
265	'''
266	# no matter the shape, transform to 3xN
NEW 267	points = np.array(points)	×
NEW 268	if points.shape[0] != 3:	×
NEW 269	points = points.T	×
270
NEW 271	x_min = np.min(points[0])	×
NEW 272	x_max = np.max(points[0])	×
NEW 273	y_min = np.min(points[1])	×
NEW 274	y_max = np.max(points[1])	×
NEW 275	z_min = np.min(points[2])	×
NEW 276	z_max = np.max(points[2])	×
277
NEW 278	return (x_min, x_max, y_min, y_max, z_min, z_max)	×
279
280	###################################################
281	# Mesh generation for different turbine components
282	###################################################
283
NEW 284	def bladeMesh(bladeData, airfoils):	×
285	'''
286	Generate mesh for blade component
287	'''
288	# Clear previous arrays
NEW 289	x = np.array([])	×
NEW 290	y = np.array([])	×
NEW 291	z = np.array([])	×
292
293	# interpolate the chord, twist, pitch axis and reference axis to the locations of the airfoils
NEW 294	chord = np.interp(bladeData['outer_shape_bem']['airfoil_position']['grid'],	×
295	bladeData['outer_shape_bem']['chord']['grid'],
296	bladeData['outer_shape_bem']['chord']['values'])
297
NEW 298	twist = np.interp(bladeData['outer_shape_bem']['airfoil_position']['grid'],	×
299	bladeData['outer_shape_bem']['twist']['grid'],
300	bladeData['outer_shape_bem']['twist']['values'])
301
NEW 302	pitch_axis = np.interp(bladeData['outer_shape_bem']['airfoil_position']['grid'],	×
303	bladeData['outer_shape_bem']['pitch_axis']['grid'],
304	bladeData['outer_shape_bem']['pitch_axis']['values'])
305
NEW 306	ref_axis_x = np.interp(bladeData['outer_shape_bem']['airfoil_position']['grid'],	×
307	bladeData['outer_shape_bem']['reference_axis']['x']['grid'],
308	bladeData['outer_shape_bem']['reference_axis']['x']['values'])
309
NEW 310	ref_axis_y = np.interp(bladeData['outer_shape_bem']['airfoil_position']['grid'],	×
311	bladeData['outer_shape_bem']['reference_axis']['y']['grid'],
312	bladeData['outer_shape_bem']['reference_axis']['y']['values'])
313
NEW 314	ref_axis_z = np.interp(bladeData['outer_shape_bem']['airfoil_position']['grid'],	×
315	bladeData['outer_shape_bem']['reference_axis']['z']['grid'],
316	bladeData['outer_shape_bem']['reference_axis']['z']['values'])
317
318
319	# extract the airfoil names at the grid locations
NEW 320	airfoil_names = [bladeData['outer_shape_bem']['airfoil_position']['labels'][i] for i in range(len(bladeData['outer_shape_bem']['airfoil_position']['grid']))]	×
321
322	# Arrays to store points
NEW 323	x = np.array([])	×
NEW 324	y = np.array([])	×
NEW 325	z = np.array([])	×
326
327	# For each blade station, we transform the coordinates and store
NEW 328	for i in range(len(ref_axis_z)):	×
329	# Get the airfoil coordinates for this section
NEW 330	af_coordinates = np.array([airfoils[airfoil_names[i]]['coordinates']['x'], airfoils[airfoil_names[i]]['coordinates']['y']]).T	×
331
332	# Scale by chord
NEW 333	af_coordinates[:, 0] = (af_coordinates[:, 0] - pitch_axis[i]) * chord[i]	×
NEW 334	af_coordinates[:, 1] = af_coordinates[:, 1] * chord[i]	×
335
336	# Create rotation matrix for twist angle
NEW 337	cos_t = np.cos(np.radians(twist[i]))	×
NEW 338	sin_t = np.sin(np.radians(twist[i]))	×
NEW 339	rot_matrix = np.array([[cos_t, -sin_t], [sin_t, cos_t]])	×
340
341	# Apply rotation
NEW 342	af_coordinates = np.dot(af_coordinates, rot_matrix.T)	×
343
344	# Translate to reference axis position
NEW 345	af_coordinates_3d = np.zeros((af_coordinates.shape[0], 3))	×
NEW 346	af_coordinates_3d[:, 0] = af_coordinates[:, 0] + ref_axis_y[i]	×
NEW 347	af_coordinates_3d[:, 1] = af_coordinates[:, 1] + -ref_axis_x[i]	×
NEW 348	af_coordinates_3d[:, 2] = np.full(af_coordinates.shape[0], ref_axis_z[i])	×
349
350	# Append points
NEW 351	x = np.append(x, af_coordinates_3d[:,0])	×
NEW 352	y = np.append(y, af_coordinates_3d[:,1])	×
NEW 353	z = np.append(z, af_coordinates_3d[:,2])	×
354
NEW 355	points = (x, y, z)	×
356
NEW 357	return points	×
358
359
NEW 360	def towerMesh(towerData):	×
361
362	# Clear previous arrays
NEW 363	x = np.array([])	×
NEW 364	y = np.array([])	×
NEW 365	z = np.array([])	×
366
NEW 367	towerGrid = np.interp(towerData['outer_shape_bem']['outer_diameter']['grid'],	×
368	towerData['outer_shape_bem']['reference_axis']['z']['grid'],
369	towerData['outer_shape_bem']['reference_axis']['z']['values'])
NEW 370	TowerOD = towerData['outer_shape_bem']['outer_diameter']['values']	×
371
372	# Create points for each cylinder segment
NEW 373	for i in range(len(towerGrid)-1):	×
374	# Get parameters for this segment
NEW 375	h_start = towerGrid[i]	×
NEW 376	h_end = towerGrid[i+1]	×
NEW 377	r_start = TowerOD[i]/2	×
NEW 378	r_end = TowerOD[i+1]/2	×
379
380	# Create points along height
NEW 381	h_points = np.linspace(h_start, h_end, 20)	×
382
383	# For each height, create a circle of points
NEW 384	for h in h_points:	×
385	# Calculate interpolated radius at this height
NEW 386	r = r_start + (r_end - r_start) * (h - h_start)/(h_end - h_start)	×
387
388	# Create circle of points at this height
NEW 389	theta = np.linspace(0, 2*np.pi, 36)	×
NEW 390	x = np.append(x, r * np.cos(theta))	×
NEW 391	y = np.append(y, r * np.sin(theta))	×
NEW 392	z = np.append(z, np.full_like(theta, h))	×
393
394
NEW 395	points = (x, y, z)	×
NEW 396	return points	×
397
398	# def hubMesh(hubData):
399
400	# # Clear previous arrays
401	# x = np.array([])
402	# y = np.array([])
403	# z = np.array([])
404
405	# dia = hubData['diameter']
406	# h = 0.5 * dia
407	# '''
408	# let the height of the hub be 0.5 of the diameter
409	# The radius of the hub along the height varies from diameter to zero at the top and follows a parabolic function
410	# '''
411	# numPoints = 100
412	# theta = np.linspace(0, 2*np.pi, numPoints)
413	# z_dists = np.linspace(0, dia/2, 20)
414
415	# for z_dist in z_dists:
416	# r = dia/2 * np.sqrt(1 - (2z_dist/dia)*2)
417	# x = np.append(x, r * np.cos(theta))
418	# y = np.append(y, r * np.sin(theta))
419	# z = np.append(z, np.full_like(theta, z_dist))
420
421	# points = (x, y, z)
422
423	# return points
424
NEW 425	def hubMesh(hubData):	×
426	# Clear previous arrays
NEW 427	x = np.array([])	×
NEW 428	y = np.array([])	×
NEW 429	z = np.array([])	×
430
NEW 431	dia = 1.25 * hubData['diameter'] # Increased diameter multiplier	×
NEW 432	total_length = dia #* 1.2 # Increased total length for better coverage	×
NEW 433	offset = -0.2 * dia # Offset to move the hub backwards	×
434
NEW 435	'''	×
436	Hub consists of a parabolic section that:
437	1. Starts slightly behind z=0 to better cover blade roots
438	2. Has maximum radius at z=0.3*length
439	3. Tapers to a point at the end
440	'''
441
NEW 442	numPoints = 100	×
NEW 443	theta = np.linspace(0, 2*np.pi, numPoints)	×
NEW 444	z_points = np.linspace(offset, total_length + offset, 30)	×
445
NEW 446	for z_dist in z_points:	×
447	# Modified parabolic profile
NEW 448	normalized_z = (z_dist - offset)/(total_length)	×
449	# Profile peaks at z=0.3*length and tapers at both ends
NEW 450	r = dia/2 * (1 - ((normalized_z - 0.3)/0.7)**2)	×
NEW 451	r = max(0, r) # Ensure radius doesn't go negative	×
452
NEW 453	x = np.append(x, r * np.cos(theta))	×
NEW 454	y = np.append(y, r * np.sin(theta))	×
NEW 455	z = np.append(z, np.full_like(theta, z_dist))	×
456
NEW 457	points = (x, y, z)	×
NEW 458	return points	×
459
460
NEW 461	def monopileMesh(monopileData):	×
462
463	# Reusing the tower mesh function for monopile
NEW 464	points = towerMesh(monopileData)	×
465
NEW 466	return points	×
467
NEW 468	def nacelleMesh(nacelleData, hubData):	×
469
470	# Clear previous arrays
NEW 471	x = np.array([])	×
NEW 472	y = np.array([])	×
NEW 473	z = np.array([])	×
474
475	# Add an offset for better visual alignment
NEW 476	xOffset = -0.2 # Shift nacelle back by to better align with hub	×
477
478	# if the nacelle dict has length, width and height, we can create a box
NEW 479	if 'length' in nacelleData['drivetrain'].keys():	×
NEW 480	length = nacelleData['drivetrain']['length']	×
NEW 481	width = nacelleData['drivetrain']['width']	×
NEW 482	height = nacelleData['drivetrain']['height']	×
NEW 483	x = np.array([length/2, length/2, -length/2, -length/2, length/2, length/2, -length/2, -length/2]).T + xOffset * length	×
NEW 484	y = np.array([width/2, -width/2, -width/2, width/2, width/2, -width/2, -width/2, width/2]).T	×
NEW 485	z = np.array([-height/2, -height/2, -height/2, -height/2, height/2, height/2, height/2, height/2]).T	×
486	else:
487	# if not, we use the overhang, distance from tower top to hub, and hub diameter to create the box such that
488	# height = 1.5 * distance from tower top to hub
489	# length = 1 * overhang
490	# width = 1 * hub diameter
NEW 491	height = 1.5 * nacelleData['drivetrain']['distance_tt_hub']	×
NEW 492	length = 1 * nacelleData['drivetrain']['overhang']	×
NEW 493	width = 1 * hubData['diameter']	×
NEW 494	x = np.array([length/2, length/2, -length/2, -length/2, length/2, length/2, -length/2, -length/2]).T + xOffset * length	×
NEW 495	y = np.array([width/2, -width/2, -width/2, width/2, width/2, -width/2, -width/2, width/2]).T	×
NEW 496	z = np.array([-height/2, -height/2, -height/2, -height/2, height/2, height/2, height/2, height/2]).T	×
497
NEW 498	points = (x, y, z)	×
499
NEW 500	return points	×
501
NEW 502	def floatingPlatformMesh(semisubData, mooringData=None):	×
503
504	# convert joints from array to dictionary
NEW 505	semisubData['joints'] = {joint['name']: joint for joint in semisubData['joints']}	×
506
507	# creating a few sub functions
NEW 508	def cylind2cart(p):	×
509	# r, theta, z
NEW 510	x = np.array([0, 0, 0])	×
NEW 511	x[0] = p[0] * np.cos(p[1])	×
NEW 512	x[1] = p[0] * np.sin(p[1])	×
NEW 513	x[2] = p[2]	×
NEW 514	return x	×
515
NEW 516	def getJointPoint(jointName):	×
517
NEW 518	joint = semisubData['joints'][jointName]	×
519
520	# search if joint has cylindrical key
NEW 521	if 'cylindrical' in joint.keys():	×
NEW 522	if joint['cylindrical']:	×
NEW 523	return cylind2cart(joint['location'])	×
524	else:
NEW 525	return joint['location']	×
526	else:
NEW 527	return joint['location']	×
528
529	# def meshCone(joint1, joint2, grid, value):
530	# grid = np.interp(grid, [0, 1], [joint1, joint2])
531
532
533	# looping through the members
NEW 534	for idx, member in enumerate(semisubData['members']):	×
535	# fetch the joints
NEW 536	joint1 = np.array(getJointPoint(member['joint1']))	×
NEW 537	joint2 = np.array(getJointPoint(member['joint2']))	×
538
539	# direction vector normalized
NEW 540	direction = (joint2 - joint1) / np.linalg.norm(joint2 - joint1)	×
541
542
NEW 543	if member['outer_shape']['shape'] == 'circular':	×
544	# create a cylinder, if last and first diameters are the same, then use pv.Cylinder
545	# else use towerMesh to create a cylinder
NEW 546	if member['outer_shape']['outer_diameter']['values'][0] == member['outer_shape']['outer_diameter']['values'][-1]:	×
NEW 547	center = joint1 + 0.5 * (joint2 - joint1) # Calculate midpoint	×
NEW 548	memberMesh = pv.Cylinder(center=center,	×
549	direction=direction,
550	height=np.linalg.norm(joint2 - joint1),
551	radius=member['outer_shape']['outer_diameter']['values'][0]/2)
552	else:
553	# Create a tapered cylinder using towerMesh-like approach
NEW 554	towerData = {	×
555	'outer_shape_bem': {
556	'outer_diameter': {
557	'grid': member['outer_shape']['outer_diameter']['grid'],
558	'values': member['outer_shape']['outer_diameter']['values']
559	},
560	'reference_axis': {
561	'z': {
562	'grid': member['outer_shape']['outer_diameter']['grid'],
563	'values': np.linspace(0, np.linalg.norm(joint2 - joint1), len(member['outer_shape']['outer_diameter']['grid']))
564	}
565	}
566	}
567	}
568
NEW 569	x, y, z = towerMesh(towerData)	×
570
571	# Rotate and translate the points to align with member direction
NEW 572	points = np.vstack((x, y, z))	×
573
574	# Calculate rotation matrix from [0,0,1] to direction vector
NEW 575	v = np.array([0, 0, 1])	×
NEW 576	rot_axis = np.cross(v, direction)	×
NEW 577	rot_angle = np.arccos(np.dot(v, direction))	×
578
NEW 579	if np.linalg.norm(rot_axis) > 0:	×
NEW 580	rot_axis = rot_axis / np.linalg.norm(rot_axis)	×
NEW 581	rotMatrix = pv.transformations.axis_angle_rotation(rot_axis, np.degrees(rot_angle))	×
NEW 582	points = np.dot(rotMatrix[:3, :3], points)	×
583
584	# Translate to joint1
NEW 585	points = points + joint1.reshape(3, 1)	×
586
NEW 587	memberMesh = pv.PolyData(points.T)	×
NEW 588	memberMesh = memberMesh.delaunay_3d()	×
589
NEW 590	elif member['outer_shape']['shape'] == 'rectangular':	×
NEW 591	'''	×
592	Definition of Rectangular member:
593	shape: rectangular
594	side_length_a:
595	grid: [0.0, 1.0]
596	values: [12.5, 12.5]
597	side_length_b:
598	grid: [0.0, 1.0]
599	values: [7.0, 7.0]
600	'''
601	# For rectangular members, create a box using the side lengths
602	# Get side lengths at grid points
NEW 603	side_a = member['outer_shape']['side_length_a']['values']	×
NEW 604	side_b = member['outer_shape']['side_length_b']['values']	×
605
606	# Calculate the member direction and perpendicular vectors
NEW 607	direction = (joint2 - joint1) / np.linalg.norm(joint2 - joint1)	×
NEW 608	perpendicular1 = np.array([-direction[1], direction[0], 0])	×
NEW 609	if np.linalg.norm(perpendicular1) < 1e-10:	×
NEW 610	perpendicular1 = np.array([1, 0, 0])	×
NEW 611	perpendicular1 = perpendicular1 / np.linalg.norm(perpendicular1)	×
NEW 612	perpendicular2 = np.cross(direction, perpendicular1)	×
613
614	# Create corners based on joint points
NEW 615	corners = np.array([	×
616	# Bottom face
617	joint1 + (-side_a[0]/2 * perpendicular1) + (-side_b[0]/2 * perpendicular2),
618	joint1 + (side_a[0]/2 * perpendicular1) + (-side_b[0]/2 * perpendicular2),
619	joint1 + (side_a[0]/2 * perpendicular1) + (side_b[0]/2 * perpendicular2),
620	joint1 + (-side_a[0]/2 * perpendicular1) + (side_b[0]/2 * perpendicular2),
621	# Top face
622	joint2 + (-side_a[1]/2 * perpendicular1) + (-side_b[1]/2 * perpendicular2),
623	joint2 + (side_a[1]/2 * perpendicular1) + (-side_b[1]/2 * perpendicular2),
624	joint2 + (side_a[1]/2 * perpendicular1) + (side_b[1]/2 * perpendicular2),
625	joint2 + (-side_a[1]/2 * perpendicular1) + (side_b[1]/2 * perpendicular2)
626	])
627
628	# Create mesh from vertices
NEW 629	faces = np.array([	×
630	[4, 0, 1, 2, 3], # bottom
631	[4, 4, 5, 6, 7], # top
632	[4, 0, 1, 5, 4], # front
633	[4, 1, 2, 6, 5], # right
634	[4, 2, 3, 7, 6], # back
635	[4, 3, 0, 4, 7] # left
636	])
NEW 637	memberMesh = pv.PolyData(corners, faces)	×
638
639
640
641
642	else:
643	# warning that the shape is not recognized
NEW 644	print(f"Shape {member['outer_shape']['shape']} not recognized")	×
645
646
647	# if the member has a key named 'axial_joints', then calculate its location in space and add it to the list of joints for future reference
NEW 648	if 'axial_joints' in member.keys():	×
NEW 649	for axial_joint in member['axial_joints']:	×
NEW 650	axial_joint_location = joint1 + axial_joint['grid'] * (joint2 - joint1)	×
NEW 651	semisubData['joints'][axial_joint['name']] = {'name': axial_joint['name'],'location': axial_joint_location, 'cylindrical': False}	×
652
653
NEW 654	if idx == 0:	×
NEW 655	floatingMesh = memberMesh	×
656	else:
NEW 657	floatingMesh = floatingMesh.merge(memberMesh)	×
658
659	# floatingMesh.plot(show_edges=True)
660
661
662	# if the semisub has a mooring system, then we need to add the mooring lines
NEW 663	if mooringData:	×
664
665	# convert nodes from array to dictionary
NEW 666	mooringData['nodes'] = {node['name']: node for node in mooringData['nodes']}	×
667
668	# Looping over the mooring lines
NEW 669	for line in mooringData['lines']:	×
670	# fetch the nodes
NEW 671	np.array(getJointPoint(member['joint1']))	×
NEW 672	node1 = np.array(getJointPoint(mooringData['nodes'][line['node1']]['joint']))	×
NEW 673	node2 = np.array(getJointPoint(mooringData['nodes'][line['node2']]['joint']))	×
674
675	# check if the distance between the nodes is equal to the unstretched length of the mooring line
NEW 676	if np.linalg.norm(node2 - node1) == line['unstretched_length']:	×
NEW 677	mooringLineMesh = pv.Line(pointa=node1, pointb=node2)	×
678	else:
679	# Create catenary mooring line using parametric equations
NEW 680	num_points = 50	×
NEW 681	t = np.linspace(0, 1, num_points)	×
682	# Calculate catenary parameter (a) based on endpoints and length
NEW 683	dx = node2[0] - node1[0]	×
NEW 684	dy = node2[1] - node1[1]	×
NEW 685	dz = node2[2] - node1[2]	×
NEW 686	L = line['unstretched_length']	×
NEW 687	a = np.sqrt(L2 - dz2) / 2 # Approximate catenary parameter	×
688
689	# Generate points along catenary
NEW 690	x = node1[0] + dx * t	×
NEW 691	y = node1[1] + dy * t	×
692	# Changed the sign before 'a' to make the catenary curve downward
NEW 693	z = node1[2] + dz * t + a * (np.cosh(t - 0.5) - np.cosh(-0.5))	×
694
NEW 695	points = np.column_stack((x, y, z))	×
NEW 696	mooringLineMesh = pv.lines_from_points(points)	×
697
NEW 698	floatingMesh = floatingMesh.merge(mooringLineMesh)	×
699
NEW 700	return floatingMesh	×
701
NEW 702	def render_our_own_delaunay(points):	×
703	'''
704	Create and fill the VTK Data Object with your own data using VTK library and pyvista high level api
705
706	Reference: https://tutorial.pyvista.org/tutorial/06_vtk/b_create_vtk.html
707	https://docs.pyvista.org/examples/00-load/create-tri-surface
708	https://docs.pyvista.org/api/core/_autosummary/pyvista.polydatafilters.reconstruct_surface#pyvista.PolyDataFilters.reconstruct_surface
709	'''
710
711	# Join the points
NEW 712	x, y, z = points	×
NEW 713	values = np.c_[x.ravel(), y.ravel(), z.ravel()] # (6400, 3) where each column is x, y, z coords	×
NEW 714	coords = numpy_to_vtk(values)	×
NEW 715	cloud = pv.PolyData(coords)	×
716	# mesh = cloud.delaunay_2d() # From point cloud, apply a 2D Delaunary filter to generate a 2d surface from a set of points on a plane.
NEW 717	mesh = cloud.delaunay_3d()	×
718
NEW 719	mesh_state = to_mesh_state(mesh)	×
720
NEW 721	return mesh_state, mesh	×
722
723
NEW 724	def main():	×
725	# fetch the turbine data from the yaml file
NEW 726	with open('../../examples/06_IEA-15-240-RWT/IEA-15-240-RWT_VolturnUS-S.yaml') as file:	×
727	# with open('../../examples/06_IEA-15-240-RWT/IEA-15-240-RWT_Monopile.yaml') as file:
NEW 728	turbine_data = yaml.load(file, Loader=yaml.FullLoader)	×
729
730	# render the turbine components
731	# mesh_state, meshTurbine, extremes_all = render_turbine(turbine_data, ['blade', 'hub', 'tower', 'nacelle', 'monopile'])
NEW 732	mesh_state, meshTurbine, extremes_all = render_turbine(turbine_data, ['blade', 'hub', 'tower', 'nacelle','floating_platform'])	×
733	# mesh_state, meshTurbine, extremes_all = render_turbine(turbine_data, ['floating_platform'])
NEW 734	meshTurbine.plot(show_edges=False)	×
735
736
737	# call main if this script is run
NEW 738	if __name__ == "__main__":	×
NEW 739	main()	×

WISDEM / WEIS / 13380864473

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous