import numpy as np
def toroidal_to_cartesien(array_in:np.ndarray):
"""
Convert toroidal coordinates to Cartesian coordinates.
This function takes an array of points in toroidal coordinates and converts them to Cartesian coordinates.
The toroidal coordinates are expected to be in the format [u, v, w, phi], where:
- u is the axial coordinate,
- v is the orthoradial coordinate,
- w is the radial coordinate,
- phi is the angle described by v w.r.t. the Cartesien coordinates.
The conversion formulas are:
- X = u
- Y = w*cos(phi) - v*sin(phi)
- Z = w*sin(phi) + v*cos(phi)
Parameters:
array_in (np.ndarray): A 2D array where each row represents a point in toroidal coordinates [u, v, w, phi].
Returns:
np.ndarray: A 2D array where each row represents a point in Cartesian coordinates [X, Y, Z].
"""
# Extract data
orthor = array_in[:, 1]
radial = array_in[:, 2]
cos_phi = np.cos(array_in[:, 3])
sin_phi = np.sin(array_in[:, 3])
array_out = np.copy(array_in)
# array_out[:, 0] = array_in[:, 0]
array_out[:, 1] = radial*cos_phi - orthor*sin_phi
array_out[:, 2] = radial*sin_phi + orthor*cos_phi
return array_out[:, :3] # Return only X Y Z
def postprocess_shell_displacement(U_s:np.ndarray, nb_modes:int, angles:np.ndarray):
"""
Postprocess the shell displacement data to extract axial, ortho-radial, and radial displacements
as a function of azimuthal angle and radial position.
The input array `U_s` is expected to be an array containing the displacement components
of the middle surface of a shell structure. The structure of `U_s` is as follows:
- For each Fourier mode (from 1 to `nb_modes`), the displacements are listed in the order:
[u^i_m, v^i_m, w^i_m, u^o_m, v^o_m, w^o_m]
- The displacement for mode 0 are listed at the end in the following order: [w^o, w^i_1, w^o_1]
Here, `i` denotes the in-plane Fourier displacement and `o` denotes the out-plane Fourier displacement,
and `m` denotes the Fourier mode. The number of modes, also called `nb_modes`, is either 3 or 6.
The `angles` array contains the azimuthal angles at which the displacements are evaluated.
The function returns a 2D array where each row corresponds to a different angle.
The first three columns are the axial, ortho-radial, and radial displacements, respectively.
"""
displacement = np.array([[0, 0, 0, phi] for phi in angles])
for i, phi in enumerate(angles):
block_list = []
for m in range(2, nb_modes+1):
block = np.array([
[np.cos(m*phi), 0, 0],
[0, np.sin(m*phi), 0],
[0, 0, np.cos(m*phi)]
])
block_list.append(block)
for m in range(2, nb_modes+1):
block = np.array([
[np.sin(m*phi), 0, 0],
[0, np.cos(m*phi), 0],
[0, 0, np.sin(m*phi)]
])
block_list.append(block)
block = np.array([
[0, 0, 0],
[0, np.sin(phi), -np.cos(phi)],
[1, np.cos(phi), np.sin(phi)],
])
block_list.append(block)
displacement[i, :-1] = np.block(block_list) @ U_s
return displacement
def postprocess_beam_displacement(U_p:np.ndarray, radius:float, angles:np.ndarray):
"""
Postprocess the beam displacement data to extract axial, ortho-radial, and radial displacements
as a function of azimuthal angle and radial position.
The input array `U_p` is expected to be an array containing the displacement components
of the section of a beam structure. The structure of `U_p` is as follows:
- The beam displacements [DX, DY, DZ] are the first three components of U_p,
- The beam rotations [RX, RY, RZ] are the last three components of U_p.
The `angles` array contains the azimuthal angles at which the displacements are evaluated.
The `radius` float is the value of the middle radius of the cylindrical shell.
The function returns a 2D array where each row corresponds to a different angle.
The first three columns are the axial, ortho-radial, and radial displacements, respectively.
"""
displacement = np.array([[0, 0, 0, phi] for phi in angles])
for i, phi in enumerate(angles):
displacement[i, 0] = U_p[0] + U_p[5]*radius*np.sin(phi) - U_p[4]*radius*np.cos(phi)
displacement[i, 1] = U_p[2]*np.sin(phi) - U_p[1]*np.cos(phi) - U_p[3]*radius
displacement[i, 2] = -(U_p[2]*np.cos(phi) + U_p[1]*np.sin(phi))
return displacement
def sum_toroidal_coords(array_list:list):
"""
Sum an arbitrary number of arrays of toroidal coordinates.
This function takes a list of arrays, each representing toroidal coordinates,
and returns a single array containing the sum of the input arrays. The input
arrays are expected to have the same shape, with the last axis representing
the angle values.
Parameters:
array_list (list): A list of arrays, where each array represents a set of toroidal coordinates.
Returns:
np.ndarray: A single array containing the sum of the input arrays, with the same shape as the individual arrays.
"""
assert len(array_list) > 0, "At least one array"
array_out = np.copy(array_list[0])
if len(array_list) > 1:
for array in array_list[1:]:
array_out[..., :-1] += array[..., :-1]
return array_out
# Define global variables
# For importing and exporting
import os
SUFIX = '6m' # '3m' or '6m'
DIRPATH = os.path.dirname(os.path.abspath(__file__))
COORDS = np.load(f"{DIRPATH}/output/coords_{SUFIX}.npy")
DEPLS = np.load(f"{DIRPATH}/output/displacements_{SUFIX}.npy")
# For post-processing
RADIUS = 1
ANGLES = np.linspace(0, 2*np.pi, 120)
# For plotting
OFFSET = 2.5
VMIN, VMAX = 0, 2
NBTICKS = 11
NBCOLORS = 2*NBTICKS-1
# Loop through each frequency
from matplotlib import colors, cm, gridspec, pyplot as plt
nb_freqs, nb_coords = np.shape(DEPLS)[:2]
nb_freqs = np.min([nb_freqs, 18])
for ii in range(nb_freqs):
fig = plt.figure(figsize=(9, 4))
gs = gridspec.GridSpec(1, 2, width_ratios=[1, 2.5], figure=fig)
ax1 = fig.add_subplot(gs[0])
ax2 = fig.add_subplot(gs[1], projection='3d')
# Loop through each coordinate
XX_all, YY_all, ZZ_all, VV_all = [], [], [], []
for jj in range(nb_coords):
# Compute the total displacement
beam_disp = postprocess_beam_displacement(U_p=DEPLS[ii, jj, :6],
radius=RADIUS, angles=ANGLES)
shell_disp = postprocess_shell_displacement(U_s=DEPLS[ii, jj, 6:],
nb_modes=int(SUFIX[0]), angles=ANGLES)
total_disp = sum_toroidal_coords([beam_disp, shell_disp])
# Define undeformed cylinder in toroidal coordinates
old_cylinder = np.array([[COORDS[jj, 0], 0, RADIUS, phi] for phi in ANGLES])
# Compute deformed cylinder in cartesien coordinates
new_cylinder = toroidal_to_cartesien(sum_toroidal_coords([old_cylinder, total_disp]))
# Collect the data for later plotting
XX_all.append(new_cylinder[:, 0])
YY_all.append(new_cylinder[:, 1])
ZZ_all.append(new_cylinder[:, 2])
VV_all.append(np.linalg.norm(new_cylinder[:, 1:], axis=1)/RADIUS)
# Convert collected data to arrays for plot
XX_all = np.array(XX_all)
YY_all = np.array(YY_all)
ZZ_all = np.array(ZZ_all)
VV_all = np.array(VV_all)
# Create the scatter plot for this frequency
ax1.scatter(YY_all, ZZ_all, c=VV_all,
norm=colors.Normalize(vmin=VMIN, vmax=VMAX),
cmap=cm.get_cmap('viridis', NBCOLORS))
ax1.grid(False)
ax1.set_xlabel('Y')
ax1.set_ylabel('Z')
ax1.set_xlim([-OFFSET, OFFSET])
ax1.set_ylim([-OFFSET, OFFSET])
ax1.set_aspect('equal')
sc = ax2.scatter(XX_all, YY_all, ZZ_all, c=VV_all,
norm=colors.Normalize(vmin=VMIN, vmax=VMAX),
cmap=cm.get_cmap('viridis', NBCOLORS))
ax2.grid(False)
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
ax2.set_zlabel('Z')
ax2.set_xlim([0, 20])
ax2.set_ylim([-OFFSET, OFFSET])
ax2.set_zlim([-OFFSET, OFFSET])
ax2.set_xticks([0, 10, 20])
cbar = fig.colorbar(sc, ax=ax2, shrink=0.8, pad=0.1)
cbar.set_ticks(np.linspace(VMIN, VMAX, NBTICKS))
cbar.set_label('Gonflement normalisé')
fig.tight_layout()
fig.savefig(f"{DIRPATH}/output/tuyau_{SUFIX}_mode_{str(ii+1)}.png")
plt.close(fig)