"""Contains functions for visualizing geometry and results.
**Contains the following classes:**
None
**Contains the following functions:**
draw: Draws a solver's Airplane(s).
animate: Animates an UnsteadyRingVortexLatticeMethodSolver's Airplane(s).
plot_results_versus_time: Plots an UnsteadyRingVortexLatticeMethodSolver's loads and
load coefficients as a function of time.
log_results: Logs a solver's load and load coefficients.
"""
from __future__ import annotations
import math
import time
import matplotlib.colors
import matplotlib.pyplot as plt
import numpy as np
import pyvista as pv
import webp
from . import (
_logging,
_parameter_validation,
_transformations,
geometry,
)
from . import operating_point as operating_point_mod
from . import (
steady_horseshoe_vortex_lattice_method,
steady_ring_vortex_lattice_method,
unsteady_ring_vortex_lattice_method,
)
_logger = _logging.get_logger("output")
# Define the color and colormaps used by the visualization functions.
_sequential_color_map = "speed"
_diverging_color_map = "delta"
_wake_vortex_color = "white"
_panel_color = "chartreuse"
_streamline_color = "orchid"
_image_surface_opacity = 0.5
_image_surface_scale = 5.0
_image_reflection_mute_factor = 0.5
_image_surface_checker_size = 25
_image_surface_color_a = np.array([40, 40, 40], dtype=np.uint8)
_image_surface_color_b = np.array([80, 80, 80], dtype=np.uint8)
_plotter_background_color = "black"
_figure_background_color = "None"
_text_color = (129, 129, 129)
_text_color_normalized: tuple[float, float, float] = (
_text_color[0] / 255,
_text_color[1] / 255,
_text_color[2] / 255,
)
_text_color_surface = (220, 220, 220)
_quality = 75.0
_window_size = [1024, 768]
# For the figure lines, use the "Prism" qualitative color map from
# carto.com/carto-colors.
_prism = [
"#5F4690",
"#1D6996",
"#38A6A5",
"#0F8554",
"#73AF48",
"#EDAD08",
"#E17C05",
"#CC503E",
"#94346E",
"#6F4070",
"#994E95",
"#666666",
]
[
_drag_color,
_side_color,
_lift_color,
_roll_color,
_pitch_color,
_yaw_color,
] = _prism[3:9]
# Set constants for the color maps, scalar bars, and text boxes.
_color_map_num_sig = 3
_bar_title_font_size = 30
_bar_label_font_size = 21
_bar_width = 0.5
_bar_position_x = 0.25
_bar_position_y = 0.05
_bar_n_labels = 2
_text_max_position = (0.85, 0.075)
_text_min_position = (0.85, 0.050)
_text_speed_position = (0.05, 0.075)
_text_font_size = 11
# Set the number of markers and the marker size for the results plots.
_num_markers = 6
_marker_size = 8
# Calculate the normalized spacing between the markers for the results plots.
_marker_spacing = 1.0 / _num_markers
[docs]
def draw(
solver: (
steady_horseshoe_vortex_lattice_method.SteadyHorseshoeVortexLatticeMethodSolver
| steady_ring_vortex_lattice_method.SteadyRingVortexLatticeMethodSolver
| unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver
),
scalar_type: str | None = None,
show_streamlines: bool | np.bool_ = False,
show_wake_vortices: bool | np.bool_ = False,
save: bool | np.bool_ = False,
testing: bool | np.bool_ = False,
) -> None:
"""Draws a solver's Airplane(s).
**Citation:**
Adapted from: vlm3.draw in AeroSandbox
Author: Peter Sharpe
Date of retrieval: 03/28/2020
:param solver: The solver whose Airplane(s) will be plotted.
:param scalar_type: Determines how to color the Panels. Setting this to None colors
the Panels uniformly. If the solver has been run, it can also be "induced drag",
"side force", or "lift", which respectively use each Panel's induced drag, side
force, and lift coefficient. The default is None.
:param show_streamlines: Set this to True to show the streamlines emanating from the
back of the Wings. If True, the solver's streamlines must have already been
calculated. Can be a bool or a numpy bool and will be converted internally to a
bool. The default is False.
:param show_wake_vortices: Set this to True to show any wake RingVortices. If True,
the solver must be an UnsteadyRingVortexLatticeMethodSolver and must have
already been run. Can be a bool or a numpy bool and will be converted internally
to a bool. The default is False.
:param save: Set this to True to save the image as a WebP. It can be a bool or a
numpy bool and will be converted internally to a bool. The default is False.
:param testing: Set this to True to close the image after one second, which is
useful for running test suites. It can be a bool or a numpy bool and will be
converted internally to a bool. The default is False.
:return: None
"""
if not isinstance(
solver,
(
steady_horseshoe_vortex_lattice_method.SteadyHorseshoeVortexLatticeMethodSolver,
steady_ring_vortex_lattice_method.SteadyRingVortexLatticeMethodSolver,
unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
),
):
raise TypeError(
"solver must be a SteadyHorseshoeVortexLatticeMethodSolver, "
"SteadyRingVortexLatticeMethodSolver, "
"or UnsteadyRingVortexLatticeMethodSolver."
)
if scalar_type is not None:
if not solver.ran:
raise RuntimeError(
"solver must have run before drawing with scalar_type not None."
)
scalar_type = _parameter_validation.str_return_str(scalar_type, "scalar_type")
if scalar_type not in ("induced drag", "side force", "lift"):
raise ValueError(
'scalar_type must be None, "induced drag", "side force", or "lift".'
)
show_streamlines = _parameter_validation.boolLike_return_bool(
show_streamlines, "show_streamlines"
)
if show_streamlines and not solver.ran:
raise RuntimeError(
"solver must have run before drawing with show_streamlines set to True."
)
if show_streamlines and len(solver.gridStreamlinePoints_GP1_CgP1) == 0:
raise RuntimeError(
"solver must have streamline points calculated before drawing with "
"show_streamlines set to True."
)
show_wake_vortices = _parameter_validation.boolLike_return_bool(
show_wake_vortices, "show_wake_vortices"
)
if show_wake_vortices and not isinstance(
solver,
unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
):
raise ValueError(
"show_wake_vortices can only be True when drawing an "
"UnsteadyRingVortexLatticeMethodSolver."
)
if show_wake_vortices and not solver.ran:
raise RuntimeError(
"solver must have run before drawing with show_wake_vortices set to True."
)
save = _parameter_validation.boolLike_return_bool(save, "save")
testing = _parameter_validation.boolLike_return_bool(testing, "testing")
# Create the Plotter and set it to use parallel projection (instead of perspective).
plotter = pv.Plotter(window_size=_window_size, lighting=None)
plotter.enable_parallel_projection() # type: ignore[call-arg]
# Get the solver's geometry and OperatingPoint.
if isinstance(
solver,
unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
):
draw_step = solver.num_steps - 1
airplanes = solver.steady_problems[draw_step].airplanes
draw_operating_point = solver.steady_problems[draw_step].operating_point
qInf__E = draw_operating_point.qInf__E
# If showing wake RingVortices, get their surfaces and plot them.
if show_wake_vortices:
wake_ring_vortex_surfaces = _get_wake_ring_vortex_surfaces(
solver, draw_step
)
plotter.add_mesh(
wake_ring_vortex_surfaces,
show_edges=True,
smooth_shading=False,
color=_wake_vortex_color,
)
else:
airplanes = solver.airplanes
draw_operating_point = solver.operating_point
qInf__E = draw_operating_point.qInf__E
# Get the Panel surfaces.
panel_surfaces = _get_panel_surfaces(airplanes)
# Plot the Panels either with scalar coloring or with a uniform color.
if scalar_type in ("induced drag", "side force", "lift"):
these_scalars = _get_scalars(airplanes, scalar_type, qInf__E)
min_scalar = round(min(these_scalars), 2)
max_scalar = round(max(these_scalars), 2)
# Choose the color map and set its limits based on if the min and max scalars
# have the same sign (sequential color map) or if they have different signs
# (diverging color map).
if np.sign(np.min(these_scalars)) == np.sign(np.max(these_scalars)):
color_map = _sequential_color_map
c_min = max(
float(np.mean(these_scalars))
- _color_map_num_sig * float(np.std(these_scalars)),
float(np.min(these_scalars)),
)
c_max = min(
float(np.mean(these_scalars))
+ _color_map_num_sig * float(np.std(these_scalars)),
float(np.max(these_scalars)),
)
else:
color_map = _diverging_color_map
c_min = -_color_map_num_sig * float(np.std(these_scalars))
c_max = _color_map_num_sig * float(np.std(these_scalars))
T_reflect = draw_operating_point.surfaceReflect_T_act_GP1_CgP1
_plot_scalars(
plotter,
these_scalars,
scalar_type,
min_scalar,
max_scalar,
color_map,
c_min,
c_max,
panel_surfaces,
text_color=_text_color_surface if T_reflect is not None else _text_color,
)
else:
plotter.add_mesh(
panel_surfaces,
show_edges=True,
color=_panel_color,
smooth_shading=False,
)
T_reflect = draw_operating_point.surfaceReflect_T_act_GP1_CgP1
image_surface_mesh = None
# If an image surface is defined, add reflected geometry. The image surface plane
# is added later, after the geometry bounds are captured.
if T_reflect is not None:
mute = _image_reflection_mute_factor
muted_edge_color = _mute_color("black", mute)
# Add reflected Panel surfaces with muted coloring.
reflected_panel_surfaces = _reflect_mesh(panel_surfaces, T_reflect)
if scalar_type in ("induced drag", "side force", "lift"):
plotter.add_mesh(
reflected_panel_surfaces,
show_edges=True,
edge_color=muted_edge_color,
cmap=_mute_colormap(color_map, mute),
clim=[c_min, c_max],
scalars=these_scalars,
smooth_shading=False,
show_scalar_bar=False,
)
else:
plotter.add_mesh(
reflected_panel_surfaces,
show_edges=True,
edge_color=muted_edge_color,
color=_mute_color(_panel_color, mute),
smooth_shading=False,
)
# Add reflected wake RingVortex surfaces if they are being shown.
if show_wake_vortices:
plotter.add_mesh(
_reflect_mesh(wake_ring_vortex_surfaces, T_reflect),
show_edges=True,
edge_color=muted_edge_color,
smooth_shading=False,
color=_mute_color(_wake_vortex_color, mute),
)
# If showing streamlines, plot them.
if show_streamlines:
# Iterate through the spanwise positions in the solver's streamline point
# ndarray.
for spanwise_position in range(solver.gridStreamlinePoints_GP1_CgP1.shape[1]):
# Get the ndarray of streamline points at this spanwise position (in the
# first Airplane's geometry axes, relative to the first Airplane's CG).
stackStreamlinePoints_GP1_CgP1 = solver.gridStreamlinePoints_GP1_CgP1[
:, spanwise_position, :
]
# Iterate through the streamline points at this spanwise position.
for point_index in range(stackStreamlinePoints_GP1_CgP1.shape[0]):
# Skip the first point because it has no previous point with which
# to make a line.
if point_index != 0:
# Get the current and last point.
point = stackStreamlinePoints_GP1_CgP1[point_index, :]
last_point = stackStreamlinePoints_GP1_CgP1[point_index - 1, :]
# Add a line to make this segment of the streamline.
plotter.add_mesh(
pv.Line(
last_point,
point,
),
show_edges=True,
color=_streamline_color,
line_width=2,
smooth_shading=False,
)
# If an image surface is defined, add the reflected streamline
# segment.
if T_reflect is not None:
reflected_point = _transformations.apply_T_to_vectors(
T_reflect,
point,
has_point=True,
)
reflected_last_point = _transformations.apply_T_to_vectors(
T_reflect,
last_point,
has_point=True,
)
plotter.add_mesh(
pv.Line(
reflected_last_point,
reflected_point,
),
show_edges=True,
color=_mute_color(_streamline_color, mute),
line_width=2,
smooth_shading=False,
)
# If an image surface is defined, save the geometry bounds (which now include
# the reflected geometry but not the image surface plane), add the image surface
# plane, then fit the camera to the saved bounds so the view is not dominated by
# the much larger image surface plane. When an image surface is present, cpos is
# not passed to show() because that would trigger an auto-fit to all actors
# (including the image surface).
if T_reflect is not None:
geometry_bounds = plotter.bounds
image_surface_result = _get_image_surface_mesh_and_texture(
draw_operating_point, geometry_bounds
)
assert image_surface_result is not None
image_surface_mesh, image_surface_texture = image_surface_result
plotter.add_mesh(
image_surface_mesh,
texture=image_surface_texture,
opacity=_image_surface_opacity,
smooth_shading=True,
)
plotter.camera.position = (-1, -1, 1)
plotter.camera.focal_point = (0, 0, 0)
plotter.camera.up = (0, 0, 1)
plotter.reset_camera(bounds=geometry_bounds) # type: ignore[call-arg]
# Set the Plotter's background color.
plotter.set_background(color=_plotter_background_color) # type: ignore[call-arg]
if not testing:
# Show the Plotter so the user can adjust the camera position and window.
# When the user closes the window, the Plotter still exists. Therefore,
# it can later be saved as an image if desired.
plotter.show(
title="Orient the view, then press any key to continue.",
cpos=(-1, -1, 1) if image_surface_mesh is None else None,
full_screen=False,
auto_close=False,
)
else:
# Show the Plotter for 1 second, then proceed automatically. This is useful
# for testing.
plotter.show(
cpos=(-1, -1, 1) if image_surface_mesh is None else None,
full_screen=False,
interactive=False,
auto_close=False,
)
time.sleep(1)
# If saving, take a screenshot, convert it to a ndarray, convert that to an Image,
# and save it as a WebP.
if save:
image = webp.Image.fromarray(
np.array(
plotter.screenshot(
filename=None,
transparent_background=True,
return_img=True,
)
)
)
webp.save_image(
img=image, file_path="Draw.webp", lossless=False, quality=_quality
)
# Close all Plotters.
pv.close_all()
[docs]
def animate(
unsteady_solver: unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
scalar_type: str | None = None,
show_wake_vortices: bool | np.bool_ = False,
save: bool | np.bool_ = False,
testing: bool | np.bool_ = False,
) -> None:
"""Animates an UnsteadyRingVortexLatticeMethodSolver's Airplane(s).
:param unsteady_solver: The UnsteadyRingVortexLatticeMethodSolver whose Airplane(s)
will be animated.
:param scalar_type: Determines how to color the Panels. Setting this to None colors
the Panels uniformly. If the solver has been run, it can also be "induced drag",
"side force", or "lift", which respectively use each Panel's induced drag, side
force, and lift coefficient. The default is None.
:param show_wake_vortices: Set this to True to show any wake RingVortices. If True,
the solver must have already been run. Can be a bool or a numpy bool and will be
converted internally to a bool. The default is False.
:param save: Set this to True to save the image as a WebP. It can be a bool or a
numpy bool and will be converted internally to a bool. The default is False.
:param testing: Set this to True to start the animation after one second, which is
useful for running test suites. It can be a bool or a numpy bool and will be
converted internally to a bool. The default is False.
:return: None
"""
if not isinstance(
unsteady_solver,
(unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,),
):
raise TypeError(
"unsteady_solver must be an UnsteadyRingVortexLatticeMethodSolver."
)
if scalar_type is not None:
if not unsteady_solver.ran:
raise RuntimeError(
"unsteady_solver must have run before animating with scalar_type not "
"None."
)
scalar_type = _parameter_validation.str_return_str(scalar_type, "scalar_type")
if scalar_type not in ("induced drag", "side force", "lift"):
raise ValueError(
'scalar_type must be None, "induced drag", "side force", or "lift".'
)
show_wake_vortices = _parameter_validation.boolLike_return_bool(
show_wake_vortices, "show_wake_vortices"
)
if show_wake_vortices and not unsteady_solver.ran:
raise RuntimeError(
"unsteady_solver must have run before animating with show_wake_vortices set"
" to True."
)
save = _parameter_validation.boolLike_return_bool(save, "save")
testing = _parameter_validation.boolLike_return_bool(testing, "testing")
first_results_step = unsteady_solver.first_results_step
# Get the solver's SteadyProblems' Airplanes. This will become a list of lists,
# with the first index being the time step and the second index identifying each
# Airplane at that time step.
step_airplanes = []
for steady_problem in unsteady_solver.steady_problems:
step_airplanes.append(steady_problem.airplanes)
# Scale down the true-speed frames per second to at most 50 fps. This is the
# maximum speed at which some programs can render WebPs.
requested_fps = 1.0 / unsteady_solver.delta_time
speed = 1.0
if requested_fps > 50.0:
speed = 50.0 / requested_fps
actual_fps = float(math.floor(requested_fps * speed))
# Create the Plotter and set it to use parallel projection (instead of perspective).
plotter = pv.Plotter(window_size=_window_size, lighting=None)
plotter.enable_parallel_projection() # type: ignore[call-arg]
# Initialize values to hold the color map choice and its limits.
c_min = 0.0
c_max = 0.0
color_map: str = ""
# Initialize variables to hold the SteadyProblems' scalars and their attributes.
all_scalars = np.empty(0, dtype=float)
min_scalar = 0.0
max_scalar = 0.0
# If coloring the Panels based on scalars, gather all the scalars across all the
# time steps and Airplanes. These will be used to set the color map limits.
if scalar_type is not None:
for step_id, airplanes in enumerate(step_airplanes):
scalars_to_add = _get_scalars(
airplanes,
scalar_type,
unsteady_solver.steady_problems[step_id].operating_point.qInf__E,
)
all_scalars = np.hstack((all_scalars, scalars_to_add))
# Choose the color map and set its limits based on if the min and max scalars
# across all time steps have the same sign (sequential color map) or if they
# have different signs (diverging color map).
if np.sign(np.min(all_scalars)) == np.sign(np.max(all_scalars)):
color_map = _sequential_color_map
c_min = max(
float(np.mean(all_scalars))
- _color_map_num_sig * float(np.std(all_scalars)),
float(np.min(all_scalars)),
)
c_max = min(
float(np.mean(all_scalars))
+ _color_map_num_sig * float(np.std(all_scalars)),
float(np.max(all_scalars)),
)
else:
color_map = _diverging_color_map
c_min = -_color_map_num_sig * float(np.std(all_scalars))
c_max = _color_map_num_sig * float(np.std(all_scalars))
min_scalar = round(min(all_scalars), 2)
max_scalar = round(max(all_scalars), 2)
# Pre-compute the image surface mesh and reflection matrix from the last time
# step's geometry so that the plane is large enough to encompass the fully
# developed wake. The mesh, texture, and reflection matrix are static and reused
# for every frame. The last step's geometry bounds (including reflected geometry
# but not the image surface plane) are also saved so the camera can be fitted to
# the geometry rather than the larger image surface.
last_step = len(step_airplanes) - 1
last_step_operating_point = unsteady_solver.steady_problems[
last_step
].operating_point
T_reflect = last_step_operating_point.surfaceReflect_T_act_GP1_CgP1
animate_text_color = _text_color_surface if T_reflect is not None else _text_color
if T_reflect is not None:
last_step_panel_surfaces = _get_panel_surfaces(step_airplanes[last_step])
reflected_last_step_panel_surfaces = _reflect_mesh(
last_step_panel_surfaces, T_reflect
)
if show_wake_vortices:
last_step_wake_surfaces = _get_wake_ring_vortex_surfaces(
unsteady_solver, last_step
)
reflected_last_step_wake_surfaces = _reflect_mesh(
last_step_wake_surfaces, T_reflect
)
combined = (
last_step_panel_surfaces.merge(last_step_wake_surfaces)
.merge(reflected_last_step_panel_surfaces)
.merge(reflected_last_step_wake_surfaces)
)
image_surface_geometry_bounds = combined.bounds
else:
combined = last_step_panel_surfaces.merge(
reflected_last_step_panel_surfaces
)
image_surface_geometry_bounds = combined.bounds
image_surface_result = _get_image_surface_mesh_and_texture(
last_step_operating_point, image_surface_geometry_bounds
)
assert image_surface_result is not None
image_surface_mesh, image_surface_texture = image_surface_result
else:
image_surface_mesh = None
image_surface_texture = None
image_surface_geometry_bounds = None
# If saving the animation, add text that displays its speed.
if save:
plotter.add_text(
text="Speed: " + str(round(100 * speed)) + "%",
position=_text_speed_position,
font_size=_text_font_size,
viewport=True,
color=animate_text_color,
)
# Get the Panel surfaces of the first time step's Airplane(s).
panel_surfaces = _get_panel_surfaces(step_airplanes[0])
# Plot the first time step's Airplanes' Panels either with scalar coloring or
# with a uniform color.
if scalar_type is not None and first_results_step == 0:
these_scalars = _get_scalars(
step_airplanes[0],
scalar_type,
unsteady_solver.steady_problems[0].operating_point.qInf__E,
)
_plot_scalars(
plotter,
these_scalars,
scalar_type,
min_scalar,
max_scalar,
color_map,
c_min,
c_max,
panel_surfaces,
text_color=_text_color_surface if T_reflect is not None else _text_color,
)
else:
plotter.add_mesh(
panel_surfaces,
show_edges=True,
color=_panel_color,
smooth_shading=False,
)
# If an image surface is defined, add reflected geometry, plot the pre-computed
# plane, set the camera direction, and fit the camera to the last time step's
# geometry bounds so the view is not dominated by the much larger image surface
# plane. When an image surface is present, cpos is not passed to show() because
# that would trigger an auto-fit to all actors (including the image surface).
if T_reflect is not None:
assert image_surface_mesh is not None
mute = _image_reflection_mute_factor
muted_edge_color = _mute_color("black", mute)
muted_panel_color = _mute_color(_panel_color, mute)
muted_wake_color = _mute_color(_wake_vortex_color, mute)
if color_map:
muted_color_map = _mute_colormap(color_map, mute)
else:
muted_color_map = None
# Add reflected Panel surfaces with muted coloring.
reflected_panel_surfaces = _reflect_mesh(panel_surfaces, T_reflect)
if scalar_type is not None and first_results_step == 0:
plotter.add_mesh(
reflected_panel_surfaces,
show_edges=True,
edge_color=muted_edge_color,
cmap=muted_color_map,
clim=[c_min, c_max],
scalars=these_scalars,
smooth_shading=False,
show_scalar_bar=False,
)
else:
plotter.add_mesh(
reflected_panel_surfaces,
show_edges=True,
edge_color=muted_edge_color,
color=muted_panel_color,
smooth_shading=False,
)
# Add the image surface plane.
plotter.add_mesh(
image_surface_mesh,
texture=image_surface_texture,
opacity=_image_surface_opacity,
smooth_shading=True,
)
plotter.camera.position = (-1, -1, 1)
plotter.camera.focal_point = (0, 0, 0)
plotter.camera.up = (0, 0, 1)
plotter.reset_camera(bounds=image_surface_geometry_bounds) # type: ignore[call-arg]
# Set the Plotter's background color.
plotter.set_background(color=_plotter_background_color) # type: ignore[call-arg]
# If not testing, show the Plotter with the first time step so the user can
# orient the view. When the user presses any key, set the title back to the
# animation title and proceed. If testing, show the Plotter with the first time
# step for 1 second, and start the animation with the current window view.
if not testing:
plotter.show(
title="Orient the view, then press any key to produce the animation.",
cpos=(-1, -1, 1) if image_surface_mesh is None else None,
full_screen=False,
auto_close=False,
)
plotter.ren_win.SetWindowName("Rendering speed not to scale.")
else:
plotter.show(
title="Rendering speed not to scale.",
cpos=(-1, -1, 1) if image_surface_mesh is None else None,
full_screen=False,
interactive=False,
auto_close=False,
)
time.sleep(1)
# Start a list to hold a WebP Image of each frame. To start, take a screenshot,
# convert it to a ndarray, and convert that to an Image.
images = [
webp.Image.fromarray(
np.array(
plotter.screenshot(
transparent_background=True,
return_img=True,
)
)
)
]
# Initialize a variable to keep track of the current time step.
current_step = 1
# Begin to iterate through the Airplane(s) from the subsequent time steps.
for airplanes in step_airplanes[1:]:
# Clear the Plotter.
plotter.clear()
# Get the Panel surfaces of this time step's Airplane(s).
panel_surfaces = _get_panel_surfaces(airplanes)
# If saving the animation, add text that displays its speed.
if save:
plotter.add_text(
text="Speed: " + str(round(100 * speed)) + "%",
position=_text_speed_position,
font_size=_text_font_size,
viewport=True,
color=animate_text_color,
)
# If showing wake RingVortices, get their surfaces and plot them.
if show_wake_vortices:
wake_ring_vortex_surfaces = _get_wake_ring_vortex_surfaces(
unsteady_solver, current_step
)
plotter.add_mesh(
wake_ring_vortex_surfaces,
show_edges=True,
smooth_shading=False,
color=_wake_vortex_color,
)
# Plot the Panels either with a uniform color or, if the current time step
# has results, with scalar coloring.
if scalar_type is not None and first_results_step <= current_step:
these_scalars = _get_scalars(
airplanes,
scalar_type,
unsteady_solver.steady_problems[current_step].operating_point.qInf__E,
)
_plot_scalars(
plotter,
these_scalars,
scalar_type,
min_scalar,
max_scalar,
color_map,
c_min,
c_max,
panel_surfaces,
text_color=(
_text_color_surface if T_reflect is not None else _text_color
),
)
else:
plotter.add_mesh(
panel_surfaces,
show_edges=True,
color=_panel_color,
smooth_shading=False,
)
# If an image surface is defined, add reflected geometry and the pre-computed
# image surface plane.
if T_reflect is not None:
assert image_surface_mesh is not None
# Add reflected Panel surfaces with muted coloring.
reflected_panel_surfaces = _reflect_mesh(panel_surfaces, T_reflect)
if scalar_type is not None and first_results_step <= current_step:
plotter.add_mesh(
reflected_panel_surfaces,
show_edges=True,
edge_color=muted_edge_color,
cmap=muted_color_map,
clim=[c_min, c_max],
scalars=these_scalars,
smooth_shading=False,
show_scalar_bar=False,
)
else:
plotter.add_mesh(
reflected_panel_surfaces,
show_edges=True,
edge_color=muted_edge_color,
color=muted_panel_color,
smooth_shading=False,
)
# Add reflected wake RingVortex surfaces if they are being shown.
if show_wake_vortices:
plotter.add_mesh(
_reflect_mesh(wake_ring_vortex_surfaces, T_reflect),
show_edges=True,
edge_color=muted_edge_color,
smooth_shading=False,
color=muted_wake_color,
)
# Add the image surface plane.
plotter.add_mesh(
image_surface_mesh,
texture=image_surface_texture,
opacity=_image_surface_opacity,
smooth_shading=True,
)
# If an image surface is present, force VTK to recalculate the scalar bar
# layout. Adding the image surface mesh with opacity causes VTK's
# UnconstrainedFontSize layout to misposition the left label (PyVista
# issue #7516).
if T_reflect is not None:
for scalar_bar_actor in plotter.scalar_bars.values():
scalar_bar_actor.Modified()
plotter.render()
# If saving, append a WebP Image of this frame to the list of Images. To do
# so, take a screenshot, convert it to a ndarray, and convert that to an Image.
if save:
images.append(
webp.Image.fromarray(
np.array(
plotter.screenshot(
filename=None,
transparent_background=True,
return_img=True,
)
)
)
)
# Increment the time step tracker.
current_step += 1
# If saving, save the list of Images as an animated WebP.
if save:
# Convert the list of WebP Images to an WebP animation.
webp.save_images(
images, "Animate.webp", fps=actual_fps, lossless=False, quality=_quality
)
# Close all the Plotters.
pv.close_all()
[docs]
def plot_results_versus_time(
unsteady_solver: unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
show: bool | np.bool_ = True,
save: bool | np.bool_ = False,
) -> None:
"""Plots an UnsteadyRingVortexLatticeMethodSolver's loads and load coefficients as a
function of time.
:param unsteady_solver: The UnsteadyRingVortexLatticeMethodSolver whose loads and
load coefficients will be plotted.
:param show: Set this to True to show the plots. It can be a bool or a numpy bool
and will be converted internally to a bool. The default is True.
:param save: Set this to True to save the plots as PNGs. It can be a bool or a numpy
bool and will be converted internally to a bool. The default is False.
:return: None
"""
if not isinstance(
unsteady_solver,
unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
):
raise TypeError(
"unsteady_solver must be an " "UnsteadyRingVortexLatticeMethodSolver."
)
show = _parameter_validation.boolLike_return_bool(show, "show")
save = _parameter_validation.boolLike_return_bool(save, "save")
if not unsteady_solver.ran:
raise RuntimeError(
"unsteady_solver must have run before plotting results versus time."
)
first_results_step = unsteady_solver.first_results_step
# Get the time step characteristics. Note that the first time step (time step
# 0), occurs at 0 seconds.
num_steps = unsteady_solver.num_steps
delta_time = unsteady_solver.delta_time
num_airplanes = unsteady_solver.num_airplanes
first_results_time_step_time = delta_time * first_results_step
final_time_step_time = delta_time * (num_steps - 1)
num_steps_to_average = num_steps - first_results_step
# Create a 1D ndarray with the time at each time step where results have been
# calculated.
times = np.linspace(
first_results_time_step_time,
final_time_step_time,
num_steps_to_average,
endpoint=True,
)
# Initialize matrices to hold the loads and load coefficients at every time step
# that has results.
forces_W = np.zeros((num_airplanes, 3, num_steps_to_average), dtype=float)
forceCoefficients_W = np.zeros(
(num_airplanes, 3, num_steps_to_average), dtype=float
)
moments_W_CgP1 = np.zeros((num_airplanes, 3, num_steps_to_average), dtype=float)
momentCoefficients_W_CgP1 = np.zeros(
(num_airplanes, 3, num_steps_to_average), dtype=float
)
# Initialize a variable to track position in the results arrays.
results_step = 0
# Iterate through the time steps and add the results to their respective matrices.
for step in range(first_results_step, num_steps):
# Get the Airplanes from the SteadyProblem at this time step.
airplanes = unsteady_solver.steady_problems[step].airplanes
# Iterate through this time step's Airplanes.
for airplane_id, airplane in enumerate(airplanes):
forces_W[airplane_id, :, results_step] = airplane.forces_W
forceCoefficients_W[airplane_id, :, results_step] = (
airplane.forceCoefficients_W
)
moments_W_CgP1[airplane_id, :, results_step] = airplane.moments_W_CgP1
momentCoefficients_W_CgP1[airplane_id, :, results_step] = (
airplane.momentCoefficients_W_CgP1
)
results_step += 1
# Iterate through the Airplane ID's to plot each Airplane's figures.
for airplane_id in range(num_airplanes):
# Initialize the four figures.
force_figure, force_axes = plt.subplots()
force_coefficients_figure, force_coefficients_axes = plt.subplots()
moment_coefficients_figure, moment_coefficients_axes = plt.subplots()
moment_figure, moment_axes = plt.subplots()
# Remove all the plots' top and right spines.
force_axes.spines.right.set_visible(False)
force_axes.spines.top.set_visible(False)
force_coefficients_axes.spines.right.set_visible(False)
force_coefficients_axes.spines.top.set_visible(False)
moment_axes.spines.right.set_visible(False)
moment_axes.spines.top.set_visible(False)
moment_coefficients_axes.spines.right.set_visible(False)
moment_coefficients_axes.spines.top.set_visible(False)
# Format all the plots' spine and label colors.
force_axes.spines.bottom.set_color(_text_color_normalized)
force_axes.spines.left.set_color(_text_color_normalized)
force_axes.xaxis.label.set_color(_text_color_normalized)
force_axes.yaxis.label.set_color(_text_color_normalized)
force_coefficients_axes.spines.bottom.set_color(_text_color_normalized)
force_coefficients_axes.spines.left.set_color(_text_color_normalized)
force_coefficients_axes.xaxis.label.set_color(_text_color_normalized)
force_coefficients_axes.yaxis.label.set_color(_text_color_normalized)
moment_coefficients_axes.spines.bottom.set_color(_text_color_normalized)
moment_coefficients_axes.spines.left.set_color(_text_color_normalized)
moment_coefficients_axes.xaxis.label.set_color(_text_color_normalized)
moment_coefficients_axes.yaxis.label.set_color(_text_color_normalized)
moment_axes.spines.bottom.set_color(_text_color_normalized)
moment_axes.spines.left.set_color(_text_color_normalized)
moment_axes.xaxis.label.set_color(_text_color_normalized)
moment_axes.yaxis.label.set_color(_text_color_normalized)
# Format all the plots' tick colors.
force_axes.tick_params(axis="x", colors=_text_color_normalized)
force_axes.tick_params(axis="y", colors=_text_color_normalized)
force_coefficients_axes.tick_params(axis="x", colors=_text_color_normalized)
force_coefficients_axes.tick_params(axis="y", colors=_text_color_normalized)
moment_coefficients_axes.tick_params(axis="x", colors=_text_color_normalized)
moment_coefficients_axes.tick_params(axis="y", colors=_text_color_normalized)
moment_axes.tick_params(axis="x", colors=_text_color_normalized)
moment_axes.tick_params(axis="y", colors=_text_color_normalized)
# Format all the plots' background colors.
force_figure.patch.set_facecolor(_figure_background_color)
force_axes.set_facecolor(_figure_background_color)
force_coefficients_figure.patch.set_facecolor(_figure_background_color)
force_coefficients_axes.set_facecolor(_figure_background_color)
moment_figure.patch.set_facecolor(_figure_background_color)
moment_axes.set_facecolor(_figure_background_color)
moment_coefficients_figure.patch.set_facecolor(_figure_background_color)
moment_coefficients_axes.set_facecolor(_figure_background_color)
# Populate the plots.
force_axes.plot(
times,
-forces_W[airplane_id, 0],
label="Induced Drag",
color=_drag_color,
marker=".",
markevery=(_marker_spacing * 0 / 3, _marker_spacing),
markersize=_marker_size,
)
force_axes.plot(
times,
forces_W[airplane_id, 1],
label="Side Force",
color=_side_color,
marker=".",
markevery=(_marker_spacing * 1 / 3, _marker_spacing),
markersize=_marker_size,
)
force_axes.plot(
times,
-forces_W[airplane_id, 2],
label="Lift",
color=_lift_color,
marker=".",
markevery=(_marker_spacing * 2 / 3, _marker_spacing),
markersize=_marker_size,
)
force_coefficients_axes.plot(
times,
-forceCoefficients_W[airplane_id, 0],
label="Induced Drag Coefficient",
color=_drag_color,
marker=".",
markevery=(_marker_spacing * 0 / 3, _marker_spacing),
markersize=_marker_size,
)
force_coefficients_axes.plot(
times,
forceCoefficients_W[airplane_id, 1],
label="Side Force Coefficient",
color=_side_color,
marker=".",
markevery=(_marker_spacing * 1 / 3, _marker_spacing),
markersize=_marker_size,
)
force_coefficients_axes.plot(
times,
-forceCoefficients_W[airplane_id, 2],
label="Lift Coefficient",
color=_lift_color,
marker=".",
markevery=(_marker_spacing * 2 / 3, _marker_spacing),
markersize=_marker_size,
)
moment_axes.plot(
times,
moments_W_CgP1[airplane_id, 0],
label="Roll",
color=_roll_color,
marker=".",
markevery=(_marker_spacing * 0 / 3, _marker_spacing),
markersize=_marker_size,
)
moment_axes.plot(
times,
moments_W_CgP1[airplane_id, 1],
label="Pitch",
color=_pitch_color,
marker=".",
markevery=(_marker_spacing * 1 / 3, _marker_spacing),
markersize=_marker_size,
)
moment_axes.plot(
times,
moments_W_CgP1[airplane_id, 2],
label="Yaw",
color=_yaw_color,
marker=".",
markevery=(_marker_spacing * 2 / 3, _marker_spacing),
markersize=_marker_size,
)
moment_coefficients_axes.plot(
times,
momentCoefficients_W_CgP1[airplane_id, 0],
label="Roll Coefficient",
color=_roll_color,
marker=".",
markevery=(_marker_spacing * 0 / 3, _marker_spacing),
markersize=_marker_size,
)
moment_coefficients_axes.plot(
times,
momentCoefficients_W_CgP1[airplane_id, 1],
label="Pitch Coefficient",
color=_pitch_color,
marker=".",
markevery=(_marker_spacing * 1 / 3, _marker_spacing),
markersize=_marker_size,
)
moment_coefficients_axes.plot(
times,
momentCoefficients_W_CgP1[airplane_id, 2],
label="Yaw Coefficient",
color=_yaw_color,
marker=".",
markevery=(_marker_spacing * 2 / 3, _marker_spacing),
markersize=_marker_size,
)
# Find and format this Airplane's name for use in the plot titles.
airplane_name = unsteady_solver.steady_problems[0].airplanes[airplane_id].name
force_title = airplane_name + " Forces vs. Time"
force_coefficient_title = airplane_name + " Force Coefficients vs. Time"
moment_title = airplane_name + " Moments vs. Time"
moment_coefficient_title = airplane_name + " Moment Coefficients vs. Time"
# Name the plots' axis labels and titles.
force_axes.set_xlabel("Time (s)", color=_text_color_normalized)
force_axes.set_ylabel("Force (N)", color=_text_color_normalized)
force_axes.set_title(force_title, color=_text_color_normalized)
force_coefficients_axes.set_xlabel("Time (s)", color=_text_color_normalized)
force_coefficients_axes.set_ylabel(
"Force Coefficient", color=_text_color_normalized
)
force_coefficients_axes.set_title(
force_coefficient_title, color=_text_color_normalized
)
moment_axes.set_xlabel("Time (s)", color=_text_color_normalized)
moment_axes.set_ylabel("Moment (N m)", color=_text_color_normalized)
moment_axes.set_title(moment_title, color=_text_color_normalized)
moment_coefficients_axes.set_xlabel("Time (s)", color=_text_color_normalized)
moment_coefficients_axes.set_ylabel(
"Moment Coefficient", color=_text_color_normalized
)
moment_coefficients_axes.set_title(
moment_coefficient_title, color=_text_color_normalized
)
# Format the plots' legends.
force_axes.legend(
facecolor=_figure_background_color,
edgecolor=_figure_background_color,
labelcolor=_text_color_normalized,
)
force_coefficients_axes.legend(
facecolor=_figure_background_color,
edgecolor=_figure_background_color,
labelcolor=_text_color_normalized,
)
moment_axes.legend(
facecolor=_figure_background_color,
edgecolor=_figure_background_color,
labelcolor=_text_color_normalized,
)
moment_coefficients_axes.legend(
facecolor=_figure_background_color,
edgecolor=_figure_background_color,
labelcolor=_text_color_normalized,
)
# Save the figures as PNGs if the user wants to do so.
if save:
force_figure.savefig(
airplane_name + " Forces.png",
dpi=300,
)
force_coefficients_figure.savefig(
airplane_name + " Force Coefficients.png",
dpi=300,
)
moment_figure.savefig(
airplane_name + " Moments.png",
dpi=300,
)
moment_coefficients_figure.savefig(
airplane_name + " Moment Coefficients.png",
dpi=300,
)
# If the user wants to show the plots, do so. This is done outside the loop so
# that plt.show() is only called once after all figures are created.
if show:
plt.show()
else:
plt.close("all")
# TEST: Consider adding integration tests for this function.
[docs]
def log_results(
solver: (
steady_horseshoe_vortex_lattice_method.SteadyHorseshoeVortexLatticeMethodSolver
| steady_ring_vortex_lattice_method.SteadyRingVortexLatticeMethodSolver
| unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver
),
) -> None:
"""Logs a solver's load and load coefficients.
The logging level must be set to INFO or lower in order to see results. See
set_up_logging for details on configuring the logging level.
:param solver: The solver whose load and load coefficients will be logged.
:return: None
"""
if isinstance(
solver,
(
steady_horseshoe_vortex_lattice_method.SteadyHorseshoeVortexLatticeMethodSolver,
steady_ring_vortex_lattice_method.SteadyRingVortexLatticeMethodSolver,
),
):
these_airplanes = solver.airplanes
solver_type = "steady"
elif isinstance(
solver,
unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
):
these_airplanes = solver.current_airplanes
if solver.unsteady_problem.movement.static:
solver_type = "static geometry unsteady"
else:
solver_type = "variable geometry unsteady"
else:
raise TypeError(
"solver must be a SteadyHorseshoeVortexLatticeMethodSolver, "
"a SteadyRingVortexLatticeMethodSolver, or an "
"UnsteadyRingVortexLatticeMethodSolver."
)
if not solver.ran:
raise RuntimeError("solver must have run before logging results.")
padding_spaces = 2
col1 = [
"FX_W",
"FY_W",
"FZ_W",
"MX_W_Cg",
"MY_W_Cg",
"MZ_W_Cg",
"cFX_W",
"cFY_W",
"cFZ_W",
"cMX_W_Cg",
"cMY_W_Cg",
"cMZ_W_Cg",
]
col1 = [label + ":" for label in col1]
col1_space = max(len(elem) for elem in col1) + padding_spaces
col3 = [
"Drag",
"Side Force",
"Lift",
"Rolling Moment",
"Pitching Moment",
"Yawing Moment",
"CDi",
"CY",
"CL",
"Cl",
"Cm",
"Cn",
]
col3 = [label + ":" for label in col3]
col3_space = max(len(elem) for elem in col3) + padding_spaces
pad = " " * padding_spaces
for airplane_num, airplane in enumerate(these_airplanes):
title1: str = ""
title2: str = ""
title3: str = ""
title4: str = ""
these_forces_W: np.ndarray = np.empty(0, dtype=float)
these_moments_W_CgP1: np.ndarray = np.empty(0, dtype=float)
these_forceCoefficients_W: np.ndarray = np.empty(0, dtype=float)
these_momentCoefficients_W_CgP1: np.ndarray = np.empty(0, dtype=float)
match solver_type:
case "steady":
title1 = f"{pad}Forces (in wind axes):"
title2 = f"{pad}Moments (in wind axes, relative to the CG):"
title3 = f"{pad}Force Coefficients (in wind axes):"
title4 = f"{pad}Moment Coefficients (in wind axes, relative to the CG):"
_forces_W = airplane.forces_W
assert _forces_W is not None
these_forces_W = _forces_W
_moments_W_CgP1 = airplane.moments_W_CgP1
assert _moments_W_CgP1 is not None
these_moments_W_CgP1 = _moments_W_CgP1
_forceCoefficients_W = airplane.forceCoefficients_W
assert _forceCoefficients_W is not None
these_forceCoefficients_W = _forceCoefficients_W
_momentCoefficients_W_CgP1 = airplane.momentCoefficients_W_CgP1
assert _momentCoefficients_W_CgP1 is not None
these_momentCoefficients_W_CgP1 = _momentCoefficients_W_CgP1
case "static geometry unsteady":
assert isinstance(
solver,
unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
)
title1 = f"{pad}Final Forces (in wind axes):"
title2 = f"{pad}Final Moments (in wind axes, relative to the CG):"
title3 = f"{pad}Final Force Coefficients (in wind axes):"
title4 = (
f"{pad}Final Moment Coefficients (in wind axes, relative to "
f"the CG):"
)
these_forces_W = solver.unsteady_problem.finalForces_W[airplane_num]
these_moments_W_CgP1 = solver.unsteady_problem.finalMoments_W_CgP1[
airplane_num
]
these_forceCoefficients_W = (
solver.unsteady_problem.finalForceCoefficients_W[airplane_num]
)
these_momentCoefficients_W_CgP1 = (
solver.unsteady_problem.finalMomentCoefficients_W_CgP1[airplane_num]
)
case "variable geometry unsteady":
assert isinstance(
solver,
unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
)
title1 = f"{pad}Final Cycle-Averaged Forces (in wind axes):"
title2 = (
f"{pad}Final Cycle-Averaged Moments (in wind axes, relative "
f"to the CG):"
)
title3 = f"{pad}Final Cycle-Averaged Force Coefficients (in wind axes):"
title4 = (
f"{pad}Final Cycle-Averaged Moment Coefficients (in wind "
f"axes, relative to the CG):"
)
these_forces_W = solver.unsteady_problem.finalMeanForces_W[airplane_num]
these_moments_W_CgP1 = solver.unsteady_problem.finalMeanMoments_W_CgP1[
airplane_num
]
these_forceCoefficients_W = (
solver.unsteady_problem.finalMeanForceCoefficients_W[airplane_num]
)
these_momentCoefficients_W_CgP1 = (
solver.unsteady_problem.finalMeanMomentCoefficients_W_CgP1[
airplane_num
]
)
case _:
raise ValueError(f"Unknown solver type: {solver_type}")
col2 = [
these_forces_W[0],
these_forces_W[1],
these_forces_W[2],
these_moments_W_CgP1[0],
these_moments_W_CgP1[1],
these_moments_W_CgP1[2],
these_forceCoefficients_W[0],
these_forceCoefficients_W[1],
these_forceCoefficients_W[2],
these_momentCoefficients_W_CgP1[0],
these_momentCoefficients_W_CgP1[1],
these_momentCoefficients_W_CgP1[2],
]
col2 = [str(np.round(val, 3)) for val in col2]
col2 = [
val + " N" if i < 3 else val + " Nm" if i < 6 else val
for i, val in enumerate(col2)
]
col2_space = max(len(elem) for elem in col2) + 2 * padding_spaces
col4 = [
-these_forces_W[0],
these_forces_W[1],
-these_forces_W[2],
these_moments_W_CgP1[0],
these_moments_W_CgP1[1],
these_moments_W_CgP1[2],
-these_forceCoefficients_W[0],
these_forceCoefficients_W[1],
-these_forceCoefficients_W[2],
these_momentCoefficients_W_CgP1[0],
these_momentCoefficients_W_CgP1[1],
these_momentCoefficients_W_CgP1[2],
]
col4 = [str(np.round(val, 3)) for val in col4]
col4 = [
val + " N" if i < 3 else val + " Nm" if i < 6 else val
for i, val in enumerate(col4)
]
_logger.info(f'Airplane "{airplane.name}":')
# Display the Reynolds number for steady solvers.
if solver_type == "steady":
assert isinstance(
solver,
(
steady_horseshoe_vortex_lattice_method.SteadyHorseshoeVortexLatticeMethodSolver,
steady_ring_vortex_lattice_method.SteadyRingVortexLatticeMethodSolver,
),
)
re = solver.reynolds_numbers[airplane_num]
_logger.info(f"{pad}Reynolds Number: {re:.2e}")
for i in range(len(col1)):
if i == 0:
_logger.info(title1)
elif i == 3:
_logger.info(title2)
elif i == 6:
_logger.info(title3)
elif i == 9:
_logger.info(title4)
s = f"{2 * pad}{col1[i]:<{col1_space}}{col2[i]:<{col2_space}}{col3[i]:<{col3_space}}{col4[i]}"
_logger.info(s)
# If the results from more Airplanes are going to be logged, log a blank
# line to separate them.
if (airplane_num + 1) < solver.num_airplanes:
_logger.info("")
def _get_panel_surfaces(
airplanes: tuple[geometry.airplane.Airplane, ...],
) -> pv.PolyData:
"""Returns a PolyData representation of the Wings' Panels' surfaces associated with
all the Airplanes in a tuple of Airplanes.
:param airplanes: The tuple of Airplanes whose Wings' Panels' surfaces will be
returned.
:return: A PolyData representation of the Airplanes' Wings' Panels' surfaces.
"""
# Initialize empty ndarrays to hold the Panels' vertices and faces.
panel_vertices = np.empty((0, 3), dtype=float)
panel_faces = np.empty(0, dtype=int)
# Initialize a variable to keep track of how many Panels have been added thus far.
panel_num = 0
# Increment through the Airplanes' Wing(s).
for airplane in airplanes:
for wing in airplane.wings:
_panels = wing.panels
assert _panels is not None
# Unravel this Wing's ndarray of Panels iterate through it.
panels = np.ravel(_panels)
for panel in panels:
# Arrange this Panel's vertices and faces into ndarrays in the
# proper form to represent PolyData surfaces.
panel_vertices_to_add = np.vstack(
(
panel.Flpp_GP1_CgP1,
panel.Frpp_GP1_CgP1,
panel.Brpp_GP1_CgP1,
panel.Blpp_GP1_CgP1,
)
)
panel_face_to_add = np.array(
[
4,
(panel_num * 4),
(panel_num * 4) + 1,
(panel_num * 4) + 2,
(panel_num * 4) + 3,
],
dtype=int,
)
# Add this Panel's vertices and faces to the ndarray of all vertices
# and faces.
panel_vertices = np.vstack((panel_vertices, panel_vertices_to_add))
panel_faces = np.hstack((panel_faces, panel_face_to_add))
# Update the number of Panels.
panel_num += 1
# Return the Panels' surfaces.
return pv.PolyData(panel_vertices, panel_faces)
def _get_image_surface_mesh_and_texture(
this_operating_point: operating_point_mod.OperatingPoint,
geometry_bounds: tuple[float, float, float, float, float, float],
) -> tuple[pv.PolyData, pv.Texture] | None:
"""Returns a PolyData plane mesh and checkerboard Texture representing the image
surface, or None if no image surface is defined.
The plane is centered at the projection of the geometry bounding box center onto the
image surface, and sized proportionally to the bounding box diagonal so that it
appears large relative to the geometry.
:param this_operating_point: The OperatingPoint that may define an image surface.
:param geometry_bounds: The (xmin, xmax, ymin, ymax, zmin, zmax) bounding box of the
geometry used to determine the plane's center and size.
:return: A tuple of (PolyData plane mesh, checkerboard Texture) representing the
image surface, or None if no image surface is defined.
"""
surface_normal = this_operating_point.surfaceNormal_GP1
surface_point = this_operating_point.surfacePoint_GP1_CgP1
if surface_normal is None or surface_point is None:
return None
# Compute the bounding box center and diagonal length.
bounds = np.array(geometry_bounds, dtype=float)
bbox_center = np.array(
[
0.5 * (bounds[0] + bounds[1]),
0.5 * (bounds[2] + bounds[3]),
0.5 * (bounds[4] + bounds[5]),
],
dtype=float,
)
bbox_diagonal = float(
np.linalg.norm(
np.array(
[
bounds[1] - bounds[0],
bounds[3] - bounds[2],
bounds[5] - bounds[4],
],
dtype=float,
)
)
)
# Project the bounding box center onto the image surface to get the plane's center.
offset = np.dot(bbox_center - surface_point, surface_normal)
plane_center = bbox_center - offset * surface_normal
# Size the plane proportionally to the bounding box diagonal.
plane_size = _image_surface_scale * bbox_diagonal
mesh = pv.Plane(
center=plane_center,
direction=surface_normal,
i_size=plane_size,
j_size=plane_size,
)
# Build a checkerboard texture image. Each cell is one pixel, so a 25 x 25
# checkerboard is a 25 x 25 x 3 RGB image.
n = _image_surface_checker_size
row = np.arange(n, dtype=int)
col = np.arange(n, dtype=int)
rr, cc = np.meshgrid(row, col, indexing="ij")
is_dark = (rr + cc) % 2 == 0
image = np.where(
is_dark[:, :, np.newaxis], _image_surface_color_a, _image_surface_color_b
)
texture = pv.numpy_to_texture(image)
return mesh, texture
def _reflect_mesh(
mesh: pv.PolyData,
T_reflect: np.ndarray,
) -> pv.PolyData:
"""Returns a copy of a PolyData mesh with its points reflected across the image
surface.
:param mesh: The PolyData mesh to reflect.
:param T_reflect: A (4,4) ndarray of floats representing the active reflection
transformation matrix (in the first Airplane's geometry axes, relative to the
first Airplane's CG).
:return: A new PolyData mesh with all points reflected across the image surface.
"""
reflected = mesh.copy()
reflected.points = _transformations.apply_T_to_vectors(
T_reflect,
mesh.points,
has_point=True,
)
return reflected
def _mute_color(
color: str | tuple[float, ...],
factor: float,
) -> tuple[float, float, float]:
"""Returns a muted version of a color by linearly interpolating it toward middle
gray.
:param color: Any color that PyVista can parse (name, hex string, RGB tuple, etc.).
:param factor: The muting factor in [0, 1]. 0 means no change, 1 means fully gray.
:return: A (R, G, B) tuple of floats in [0, 1].
"""
rgb = np.array(pv.Color(color).float_rgb)
gray = np.full(3, 0.5)
muted = rgb + factor * (gray - rgb)
return float(muted[0]), float(muted[1]), float(muted[2])
def _mute_colormap(
cmap_name: str,
factor: float,
) -> matplotlib.colors.ListedColormap:
"""Returns a muted version of a named colormap by linearly interpolating each color
toward middle gray.
:param cmap_name: The name of a Matplotlib or cmocean colormap.
:param factor: The muting factor in [0, 1]. 0 means no change, 1 means fully gray.
:return: A ListedColormap with muted colors.
"""
try:
cmap = plt.get_cmap(cmap_name)
except ValueError:
cmap = plt.get_cmap("cmo." + cmap_name)
colors = cmap(np.linspace(0, 1, 256))
gray = 0.5
colors[:, :3] = colors[:, :3] + factor * (gray - colors[:, :3])
return matplotlib.colors.ListedColormap(colors)
def _get_wake_ring_vortex_surfaces(
solver: unsteady_ring_vortex_lattice_method.UnsteadyRingVortexLatticeMethodSolver,
step: int,
) -> pv.PolyData:
"""Returns the PolyData representation of the surfaces of an
UnsteadyRingVortexLatticeMethodSolver's wake RingVortices at a given time step.
:param solver: The UnsteadyRingVortexLatticeMethodSolver with the wake RingVortices
to process.
:param step: The time step at which to process the wake RingVortices.
:return: The PolyData representation of the wake RingVortices.
"""
num_wake_ring_vortices = solver.list_num_wake_vortices[step]
stackFrwrvp_GP1_CgP1 = solver.listStackFrwrvp_GP1_CgP1[step]
stackFlwrvp_GP1_CgP1 = solver.listStackFlwrvp_GP1_CgP1[step]
stackBlwrvp_GP1_CgP1 = solver.listStackBlwrvp_GP1_CgP1[step]
stackBrwrvp_GP1_CgP1 = solver.listStackBrwrvp_GP1_CgP1[step]
# Initialize empty ndarrays to hold each wake RingVortex's vertices and face.
wake_ring_vortex_vertices = np.zeros((0, 3), dtype=float)
wake_ring_vortex_faces = np.zeros(0, dtype=int)
for wake_ring_vortex_num in range(num_wake_ring_vortices):
Frwrvp_GP1_CgP1 = stackFrwrvp_GP1_CgP1[wake_ring_vortex_num]
Flwrvp_GP1_CgP1 = stackFlwrvp_GP1_CgP1[wake_ring_vortex_num]
Blwrvp_GP1_CgP1 = stackBlwrvp_GP1_CgP1[wake_ring_vortex_num]
Brwrvp_GP1_CgP1 = stackBrwrvp_GP1_CgP1[wake_ring_vortex_num]
wake_ring_vortex_vertices_to_add = np.vstack(
(
Flwrvp_GP1_CgP1,
Frwrvp_GP1_CgP1,
Brwrvp_GP1_CgP1,
Blwrvp_GP1_CgP1,
)
)
wake_ring_vortex_face_to_add = np.array(
[
4,
(wake_ring_vortex_num * 4),
(wake_ring_vortex_num * 4) + 1,
(wake_ring_vortex_num * 4) + 2,
(wake_ring_vortex_num * 4) + 3,
],
dtype=int,
)
# Stack this wake RingVortex's vertices and faces to the ndarrays of all wake
# RingVortices' vertices and faces.
wake_ring_vortex_vertices = np.vstack(
(wake_ring_vortex_vertices, wake_ring_vortex_vertices_to_add)
)
wake_ring_vortex_faces = np.hstack(
(wake_ring_vortex_faces, wake_ring_vortex_face_to_add)
)
# Increment the wake RingVortex counter.
wake_ring_vortex_num += 1
# Return the wake RingVortex surfaces.
return pv.PolyData(wake_ring_vortex_vertices, wake_ring_vortex_faces)
def _get_scalars(
airplanes: tuple[geometry.airplane.Airplane, ...],
scalar_type: str,
qInf__E: float,
) -> np.ndarray:
"""Returns the load coefficient values from a SteadyProblem's Airplanes' Wings'
Panels.
:param airplanes: The tuple of Airplanes with the scalars to return.
:param scalar_type: Determines which load coefficient to return as scalars. Can be
"induced drag", "side force", or "lift", which respectively use each Panel's
induced drag, side force, and lift coefficient.
:param qInf__E: The current freestream dynamic pressure experienced by this
SteadyProblem's Airplane(s) (observed in the Earth frame). The units are in
Pascals.
:return: A (N,) ndarray of floats representing the N Panels' load coefficients.
"""
scalars = np.empty(0, dtype=float)
# Iterate through the Airplanes' Wings.
for airplane in airplanes:
for wing in airplane.wings:
_panels = wing.panels
assert _panels is not None
# Unravel this Wing's ndarray of Panels iterate through them.
these_panels = np.ravel(_panels)
for this_panel in these_panels:
# Stack this Panel's scalars.
if scalar_type == "induced drag":
this_induced_drag_coefficient = (
-this_panel.forces_W[0] / qInf__E / this_panel.area
)
scalars = np.hstack((scalars, this_induced_drag_coefficient))
if scalar_type == "side force":
this_side_force_coefficient = (
this_panel.forces_W[1] / qInf__E / this_panel.area
)
scalars = np.hstack((scalars, this_side_force_coefficient))
if scalar_type == "lift":
this_lift_coefficient = (
-this_panel.forces_W[2] / qInf__E / this_panel.area
)
scalars = np.hstack((scalars, this_lift_coefficient))
# Return the resulting ndarray of scalars.
return scalars
def _plot_scalars(
plotter: pv.Plotter,
these_scalars: np.ndarray,
scalar_type: str,
min_scalar: float,
max_scalar: float,
color_map: str,
c_min: float,
c_max: float,
panel_surfaces: pv.PolyData,
text_color: tuple[int, int, int] = _text_color,
) -> None:
"""Plots a scalar bar, the surfaces of a set of Panels with particular scalars, and
labels for the minimum and maximum scalar values.
:param plotter: The Plotter used for visualization.
:param these_scalars: A (N,) ndarray of floats representing the N Panels' load
coefficients.
:param scalar_type: Which load coefficient is represented by the scalars. Can be
"induced drag", "side force", or "lift".
:param min_scalar: Minimum scalar value, which is displayed as text on the Plotter.
:param max_scalar: Maximum scalar value, which is displayed as text on the Plotter.
:param color_map: Name of the color map to use for scalar visualization. Check the
pyvista.add_mesh documentation for the list of acceptable values.
:param c_min: Lower bound for the color map scaling.
:param c_max: Upper bound for the color map scaling.
:param panel_surfaces: PolyData representing the Panels' surfaces.
:param text_color: The color used for the scalar bar and label text. The default is
_text_color.
:return: None
"""
scalar_bar_args = dict(
title=scalar_type.title() + " Coefficient",
title_font_size=_bar_title_font_size,
label_font_size=_bar_label_font_size,
width=_bar_width,
position_x=_bar_position_x,
position_y=_bar_position_y,
n_labels=_bar_n_labels,
fmt="%.2f",
color=text_color,
)
plotter.add_mesh(
panel_surfaces,
show_edges=True,
cmap=color_map, # type: ignore[arg-type]
clim=[c_min, c_max],
scalars=these_scalars,
smooth_shading=False,
scalar_bar_args=scalar_bar_args, # type: ignore[arg-type]
)
plotter.add_text(
text="Max: " + str(max_scalar),
position=_text_max_position,
font_size=_text_font_size,
viewport=True,
color=text_color,
)
plotter.add_text(
text="Min: " + str(min_scalar),
position=_text_min_position,
font_size=_text_font_size,
viewport=True,
color=text_color,
)