From 478145dca2b714ec57cbb5f85fc1826b455b228b Mon Sep 17 00:00:00 2001 From: Alexandre Hazo Date: Mon, 5 May 2025 13:45:26 +0200 Subject: [PATCH 1/3] Implement base transformation from local to global for stress and strain tensors (shell elements) --- vortex_radioss/animtod3plot/Anim_to_D3plot.py | 238 ++++++++++++++++-- 1 file changed, 211 insertions(+), 27 deletions(-) diff --git a/vortex_radioss/animtod3plot/Anim_to_D3plot.py b/vortex_radioss/animtod3plot/Anim_to_D3plot.py index 1da9447..03a2491 100644 --- a/vortex_radioss/animtod3plot/Anim_to_D3plot.py +++ b/vortex_radioss/animtod3plot/Anim_to_D3plot.py @@ -54,7 +54,7 @@ def part_shell_mass(*data): n4 = node_coordinates[shell_node_indexes][:,3] ux = n1[:,0] - n3[:,0]; uy = n1[:,1] - n3[:,1]; uz = n1[:,2] - n3[:,2] vx = n4[:,0] - n2[:,0]; vy = n4[:,1] - n2[:,1]; vz = n4[:,2] - n2[:,2] - i = uy*vz - uz*vy; j = uz*vx - ux*vz; k = ux*vy - uy*vx; + i = uy*vz - uz*vy; j = uz*vx - ux*vz; k = ux*vy - uy*vx area = np.sqrt((i*i)+(j*j)+(k*k))*0.5 mass = area * shell_thicknesses * shell_densities * shell_is_alive _ = np.bincount(shell_part_indexes, mass, minlength = n_parts) @@ -76,6 +76,131 @@ def sum2Scalars(*data): print(data2) return sum(data1) + sum(data2) + @staticmethod + def normalize(v): + norm = np.linalg.norm(v) + return v / norm if norm != 0 else v + + def get_radioss_local_axes(nodes): + + """ + Computes the local Radioss coordinate system (x, y, z) and express it in the global system, + based on the node coordinates of a shell element. + + Parameters: + - nodes: np.ndarray of shape (4, 3) or (3, 3), representing the coordinates + of the element nodes [n1, n2, n3, n4] or [n1, n2, n3] + + Returns: + - x, y, z: np.ndarray (3,), local orthonormal axes expressed in the global frame. + """ + + # For 4-node shell element + if nodes.shape[0] == 4: + n1, n2, n3, n4 = nodes + + # Midpoints + m14 = 0.5 * (n1 + n4) + m23 = 0.5 * (n2 + n3) + m12 = 0.5 * (n1 + n2) + m34 = 0.5 * (n3 + n4) + + # Natural system (isoparametric frame) vectors. + xi = m23 - m14 + eta = m34 - m12 + + # Normalize xi and eta to compute the angle between them + xi_n = convert.normalize(xi) + eta_n = convert.normalize(eta) + cos_alpha = np.dot(xi_n, eta_n) + + alpha = np.arccos(cos_alpha) + + # local z-axis (normal to the shell surface) + z = np.cross(xi, eta) + + # Compute x and y so that they are orthogonal, by symmetrically rotating xi and eta + # to form a 90° angle between them, while preserving the same angle between xi and x as between eta and y + + if alpha == np.pi / 2: + x = xi_n + y = eta_n + else: + theta = (np.pi / 2 - alpha) / 2 + x = np.cos(theta) * xi_n - np.sin(theta) * eta_n + y = np.cos(theta) * eta_n - np.sin(theta) * xi_n + + + # For 3-node (triangular) shell elements + elif nodes.shape[0] == 3: + n1, n2, n3 = nodes + + x = n2 - n1 + z = np.cross(n2 - n1, n3 - n1) + y = np.cross(z, x) + + # Normalize the resulting local basis vectors + x = convert.normalize(x) + y = convert.normalize(y) + z = convert.normalize(z) + + return np.array([x, y, z]) + + + @staticmethod + def get_rotation_matrix(local_axes_radioss): + """ + Convert local axes (given as row vectors) into rotation matrices for tensor transformation. + + Parameters: + - local_axes_radioss: array (N, 3, 3), local frames with vectors as rows + + Returns: + - rotation_matrices: array (N, 3, 3), local frames with vectors as columns, for R · S_local · Rᵀ + """ + # Transpose each frame to switch from row-wise to column-wise vectors + + return np.transpose(local_axes_radioss, axes=(0, 2, 1)) + + + @staticmethod + def rotate_tensor(vecs , R): + """ + Applies rotation to a batch of 2D stress or strain tensors, + considering only in-plane components (out-of-plane components are zero). + + Parameters: + - vecs: (n, 3) array of [ox, oy, oxy] in local coordinates. + - R: (n, 3, 3) array of rotation matrices (local to global). + + Returns: + - (n, 6) array of [ox, oy, oz, oxy, oyz, oxz] in global coordinates, + where oz, oyz, oxz are zero. + """ + n = vecs.shape[0] + + # Build local 3x3 tensors (n, 3, 3) for stress/strain in-plane + S_local = np.zeros((n, 3, 3)) + S_local[:, 0, 0] = vecs[:, 0] # ox + S_local[:, 1, 1] = vecs[:, 1] # oy + S_local[:, 0, 1] = S_local[:, 1, 0] = vecs[:, 2] # σxy + + # Apply rotation: S_global = R @ S_local @ R.T + Rt = np.transpose(R, axes=(0, 2, 1)) + S_global = np.matmul(np.matmul(R, S_local), Rt) + + # Extract the components and set the out-of-plane components to zero + result = np.stack([ + S_global[:, 0, 0], # ox + S_global[:, 1, 1], # oy + np.zeros(n), # oz = 0 (out-of-plane) + S_global[:, 0, 1], # oxy + np.zeros(n), # oyz = 0 (out-of-plane) + np.zeros(n) # oxz = 0 (out-of-plane) + ], axis=-1) + + return result + @staticmethod def element_shell_internal_energy(*data): @@ -89,29 +214,46 @@ def element_shell_internal_energy(*data): @staticmethod def element_shell_stress(*data): - - # Mid-surface stresses if present are not converted - # Only Upper and Lower stresses are converted - # Out of plane stresses not converted - # Co-ordinate systems not corrected for - should be OK for Von-Mises, Tresca etc and \ - # stress principles - # [σx, σy, σz, σxy, σyz, σxz] - - shell_num = len(data[0]) - nip = data[2][0] - - out = np.zeros(shape=(shell_num, nip, 6)) - - # Top integration point - out[:, -1, 0] = data[0][:, 0] - out[:, -1, 1] = data[0][:, 1] - out[:, -1, 3] = data[0][:, 2] - - # Bottom integration point - out[:, 0, 0] = data[1][:, 0] - out[:, 0, 1] = data[1][:, 1] - out[:, 0, 3] = data[1][:, 2] - + + """ + Converts local shell stresses [ox, oy, oxy] into global stresses + [ox, oy, oz, oxy, oyz, oxz] using rotation matrices and rotate_tensor. + + Mid surface stresses if present are not converted + Only Upper and Lower stresses are converted + out of plane stresses arre not converted + """ + + shell_num = len(data[0]) + nip, rotation_matrices = data[2] + + out = np.zeros((shell_num, nip, 6)) + + top = convert.rotate_tensor(np.array(data[0]), rotation_matrices) + bottom = convert.rotate_tensor(np.array(data[1]), rotation_matrices) + + out[:, -1] = top # Top integration point + out[:, 0] = bottom # Bottom integration point + + return out + + + @staticmethod + def element_shell_strain(*data): + + #same as element shell stress + + shell_num = len(data[0]) + nip, rotation_matrices = data[2] + + out = np.zeros((shell_num, nip, 6)) + + top = convert.rotate_tensor(np.array(data[0]), rotation_matrices) + bottom = convert.rotate_tensor(np.array(data[1]), rotation_matrices) + + out[:, -1] = top # Top point + out[:, 0] = bottom # Bottom point + return out @staticmethod @@ -557,11 +699,33 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent, database_extent_binary[flag][0] = [] database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_is_alive] - database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_stress] + database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_stress] + database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_strain] database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_effective_plastic_strain] database_extent_binary[flag][1] = _[0] + [ArrayType.element_shell_thickness] - database_extent_binary[flag][1] = _[1] + [ArrayType.element_shell_internal_energy] + database_extent_binary[flag][1] = _[1] + [ArrayType.element_shell_internal_energy] + + + node_coordinates = rr.arrays["node_coordinates"] # shape: (n_nodes, 3) + shell_node_indexes = rr.arrays["element_shell_node_indexes"] + + # Extract the node matrix per element (shape: n_elements x 4 x 3) + nodes_per_elem = node_coordinates[shell_node_indexes] + + # For triangular elements, the 4th node is automatically duplicated (same as the 3rd) to maintain a (4, 3) shape. + # We remove it here so that get_radioss_local_axes correctly handle triangles as 3-node elements. + + is_triangle = np.all(nodes_per_elem[:, 2, :] == nodes_per_elem[:, 3, :], axis=1) + nodes_per_elem = [nodes[:3] if is_t else nodes for nodes, is_t in zip(nodes_per_elem, is_triangle)] + + #matrix containing the local x, y, z axes of each element + local_axes_radioss = np.array([ + convert.get_radioss_local_axes(nodes) for nodes in nodes_per_elem + ]) + + # the rotation matrices from local to global coordinate systems for each element + rotation_matrices = convert.get_rotation_matrix(local_axes_radioss) # Dyna output array_requirements[ArrayType.element_shell_thickness] = {} @@ -601,8 +765,28 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent, _["shape"] = (1,n_shell, nip_shell, 6) _["convert"] = convert.element_shell_stress _["tracker"] = shell_ids_tracker - _["additional"] = [nip_shell] + _["additional"] = [nip_shell,rotation_matrices] + + # Dyna output + array_requirements[ArrayType.element_shell_strain] = {} + _ = array_requirements[ArrayType.element_shell_strain] + # Radioss outputs needed to comptute Dyna output + _["dependents"] = ["Strain (upper)","Strain (lower)"] + _["shape"] = (1,n_shell, nip_shell, 6) + _["convert"] = convert.element_shell_strain + _["tracker"] = shell_ids_tracker + _["additional"] = [nip_shell,rotation_matrices] + # Dyna output + array_requirements[ArrayType.element_shell_strain] = {} + _ = array_requirements[ArrayType.element_shell_strain] + # Radioss outputs needed to comptute Dyna output + _["dependents"] = ["Strain (upper)","Strain (lower)"] + _["shape"] = (1,n_shell, nip_shell, 6) + _["convert"] = convert.element_shell_strain + _["tracker"] = shell_ids_tracker + _["additional"] = [nip_shell,rotation_matrices] + # Dyna output array_requirements[ArrayType.element_shell_effective_plastic_strain] = {} _ = array_requirements[ArrayType.element_shell_effective_plastic_strain] From 44368a8ba893d2126e1e354ea933a904b0d3cac9 Mon Sep 17 00:00:00 2001 From: Alexandre Hazo Date: Fri, 9 May 2025 13:57:42 +0200 Subject: [PATCH 2/3] Vectorization of local axes computation for shell elements --- vortex_radioss/animtod3plot/Anim_to_D3plot.py | 256 +++++++++--------- 1 file changed, 134 insertions(+), 122 deletions(-) diff --git a/vortex_radioss/animtod3plot/Anim_to_D3plot.py b/vortex_radioss/animtod3plot/Anim_to_D3plot.py index 03a2491..76ed9d0 100644 --- a/vortex_radioss/animtod3plot/Anim_to_D3plot.py +++ b/vortex_radioss/animtod3plot/Anim_to_D3plot.py @@ -75,93 +75,11 @@ def sum2Scalars(*data): print(data1) print(data2) return sum(data1) + sum(data2) - - @staticmethod - def normalize(v): - norm = np.linalg.norm(v) - return v / norm if norm != 0 else v - - def get_radioss_local_axes(nodes): - - """ - Computes the local Radioss coordinate system (x, y, z) and express it in the global system, - based on the node coordinates of a shell element. - - Parameters: - - nodes: np.ndarray of shape (4, 3) or (3, 3), representing the coordinates - of the element nodes [n1, n2, n3, n4] or [n1, n2, n3] - - Returns: - - x, y, z: np.ndarray (3,), local orthonormal axes expressed in the global frame. - """ - - # For 4-node shell element - if nodes.shape[0] == 4: - n1, n2, n3, n4 = nodes - - # Midpoints - m14 = 0.5 * (n1 + n4) - m23 = 0.5 * (n2 + n3) - m12 = 0.5 * (n1 + n2) - m34 = 0.5 * (n3 + n4) - - # Natural system (isoparametric frame) vectors. - xi = m23 - m14 - eta = m34 - m12 - - # Normalize xi and eta to compute the angle between them - xi_n = convert.normalize(xi) - eta_n = convert.normalize(eta) - cos_alpha = np.dot(xi_n, eta_n) - - alpha = np.arccos(cos_alpha) - - # local z-axis (normal to the shell surface) - z = np.cross(xi, eta) - - # Compute x and y so that they are orthogonal, by symmetrically rotating xi and eta - # to form a 90° angle between them, while preserving the same angle between xi and x as between eta and y - - if alpha == np.pi / 2: - x = xi_n - y = eta_n - else: - theta = (np.pi / 2 - alpha) / 2 - x = np.cos(theta) * xi_n - np.sin(theta) * eta_n - y = np.cos(theta) * eta_n - np.sin(theta) * xi_n - - - # For 3-node (triangular) shell elements - elif nodes.shape[0] == 3: - n1, n2, n3 = nodes - - x = n2 - n1 - z = np.cross(n2 - n1, n3 - n1) - y = np.cross(z, x) - - # Normalize the resulting local basis vectors - x = convert.normalize(x) - y = convert.normalize(y) - z = convert.normalize(z) - - return np.array([x, y, z]) - @staticmethod - def get_rotation_matrix(local_axes_radioss): - """ - Convert local axes (given as row vectors) into rotation matrices for tensor transformation. - - Parameters: - - local_axes_radioss: array (N, 3, 3), local frames with vectors as rows - - Returns: - - rotation_matrices: array (N, 3, 3), local frames with vectors as columns, for R · S_local · Rᵀ - """ - # Transpose each frame to switch from row-wise to column-wise vectors - - return np.transpose(local_axes_radioss, axes=(0, 2, 1)) - + def normalize(v): + norm = np.linalg.norm(v, axis=1, keepdims=True) + return np.divide(v, norm, where=norm != 0, out=np.zeros_like(v)) @staticmethod def rotate_tensor(vecs , R): @@ -274,8 +192,27 @@ def element_shell_effective_plastic_strain(*data): # Bottom integration point out[:, 0] = data[1] - return out - + return out + + @staticmethod + def element_solid_stress(*data): + # Input strain data has shape (n_solids, 6), but D3plot expects (n_solids, nip_solid, 6) + solid_num = len(data[0]) + nip = data[1][0] #Number of integration points (nip_solid) + out = np.zeros((solid_num, nip, 6)) + out[:, 0, :] = data[0] + return out + + @staticmethod + def element_solid_strain(*data): + #same as element solid stress + + solid_num = len(data[0]) + nip = data[1][0] + out = np.zeros((solid_num, nip, 6)) + out[:, 0, :] = data[0] + return out + class readAndConvert: def __init__( @@ -400,9 +337,16 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent, "We sort the ids into ascending order and generate a mapping function for the scalar and tensor data" n_nodes = rr.raw_header["nbNodes"] - - n_shell = 0 + n_shell = rr.raw_header["nbFacets"] + n_solids = rr.raw_header["nbEFunc3D"] + nip_shell = 2 + + nip_solid = 1 # Number of integration points per solid element (NIP). + # This value should be either 1 or 8, depending on the model configuration. + # However, when the model uses 8 integration points, the solver currently + # provides the average value instead of individual data points. + allowable_part_strings = ["_rigid_wall_", ": RIGIDWALL_"] if rr.raw_header["nbFacets"] > 0: @@ -688,10 +632,10 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent, _["tracker"] = node_ids_tracker _["additional"] = [] - if rr.raw_header["nbEFunc"] > 0: + if rr.raw_header["nbFacets"] > 0: flag = "SHELLS" - + database_extent_binary[flag] = {} _ = database_extent_binary[flag] @@ -700,13 +644,16 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent, database_extent_binary[flag][0] = [] database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_is_alive] database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_stress] - database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_strain] database_extent_binary[flag][0] = _[0] + [ArrayType.element_shell_effective_plastic_strain] database_extent_binary[flag][1] = _[0] + [ArrayType.element_shell_thickness] database_extent_binary[flag][1] = _[1] + [ArrayType.element_shell_internal_energy] + + # WEE need to Computes the local Radioss coordinate system (x, y, z) and express it in the global system, + # based on the node coordinates of a shell element. + node_coordinates = rr.arrays["node_coordinates"] # shape: (n_nodes, 3) shell_node_indexes = rr.arrays["element_shell_node_indexes"] @@ -717,15 +664,65 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent, # We remove it here so that get_radioss_local_axes correctly handle triangles as 3-node elements. is_triangle = np.all(nodes_per_elem[:, 2, :] == nodes_per_elem[:, 3, :], axis=1) - nodes_per_elem = [nodes[:3] if is_t else nodes for nodes, is_t in zip(nodes_per_elem, is_triangle)] - #matrix containing the local x, y, z axes of each element - local_axes_radioss = np.array([ - convert.get_radioss_local_axes(nodes) for nodes in nodes_per_elem - ]) - - # the rotation matrices from local to global coordinate systems for each element - rotation_matrices = convert.get_rotation_matrix(local_axes_radioss) + # Separate the triangles and quadrilaterals + triangles = nodes_per_elem[is_triangle, :3] # shape (n_tri, 3, 3) + quads = nodes_per_elem[~is_triangle] # shape (n_quad, 4, 3) + + local_axes = np.zeros((nodes_per_elem.shape[0], 3, 3)) # (n_elem, 3, 3) + + # --- Vectorization for TRIANGLES --- + if triangles.shape[0] > 0: + n1, n2, n3 = triangles[:, 0], triangles[:, 1], triangles[:, 2] + x = n2 - n1 + z = np.cross(n2 - n1, n3 - n1) + y = np.cross(z, x) + + # Normalize the resulting local basis vectors + x = convert.normalize(x) + y = convert.normalize(y) + z = convert.normalize(z) + + local_axes[is_triangle] = np.stack([x, y, z], axis=1) + + # --- Vectorization for 4 Nodes elements --- + if quads.shape[0] > 0: + n1, n2, n3, n4 = quads[:, 0], quads[:, 1], quads[:, 2], quads[:, 3] + + m14 = 0.5 * (n1 + n4) + m23 = 0.5 * (n2 + n3) + m12 = 0.5 * (n1 + n2) + m34 = 0.5 * (n3 + n4) + + # Natural system (isoparametric frame) vectors. + xi = m23 - m14 + eta = m34 - m12 + + # Normalize xi and eta to compute the angle between them + xi_n = convert.normalize(xi) + eta_n = convert.normalize(eta) + cos_alpha = np.einsum('ij,ij->i', xi_n, eta_n) + alpha = np.arccos(np.clip(cos_alpha, -1.0, 1.0)) + + # local z-axis (normal to the shell surface) + z = np.cross(xi, eta) + z = convert.normalize(z) + + theta = (np.pi / 2 - alpha) / 2 + theta = theta[:, None] + + # Compute x and y so that they are orthogonal, by symmetrically rotating xi and eta + # to form a 90° angle between them, while preserving the same angle between xi and x as between eta and y + x = np.cos(theta) * xi_n - np.sin(theta) * eta_n + y = np.cos(theta) * eta_n - np.sin(theta) * xi_n + + # Normaliser x et y + x = convert.normalize(x) + y = convert.normalize(y) + + local_axes[~is_triangle] = np.stack([x, y, z], axis=1) + + rotation_matrices = np.transpose(local_axes, axes=(0, 2, 1)) # Dyna output array_requirements[ArrayType.element_shell_thickness] = {} @@ -767,16 +764,27 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent, _["tracker"] = shell_ids_tracker _["additional"] = [nip_shell,rotation_matrices] - # Dyna output - array_requirements[ArrayType.element_shell_strain] = {} - _ = array_requirements[ArrayType.element_shell_strain] + # Dyna output + array_requirements[ArrayType.element_shell_effective_plastic_strain] = {} + _ = array_requirements[ArrayType.element_shell_effective_plastic_strain] # Radioss outputs needed to comptute Dyna output - _["dependents"] = ["Strain (upper)","Strain (lower)"] - _["shape"] = (1,n_shell, nip_shell, 6) - _["convert"] = convert.element_shell_strain - _["tracker"] = shell_ids_tracker - _["additional"] = [nip_shell,rotation_matrices] - + _["dependents"] = ["element_shell_plastic_strain_upper", 'element_shell_plastic_strain_lower'] + _["shape"] = (1, n_shell, nip_shell) + _["convert"] = convert.element_shell_effective_plastic_strain + _["tracker"] = shell_ids_tracker + _["additional"] = [nip_shell] + + + flag = "STRFLG" + + database_extent_binary[flag] = {} + _ = database_extent_binary[flag] + + database_extent_binary[flag][0] = [] + + if n_shell > 0: + database_extent_binary[flag][0] += [ArrayType.element_shell_strain] + # Dyna output array_requirements[ArrayType.element_shell_strain] = {} _ = array_requirements[ArrayType.element_shell_strain] @@ -787,23 +795,27 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent, _["tracker"] = shell_ids_tracker _["additional"] = [nip_shell,rotation_matrices] - # Dyna output - array_requirements[ArrayType.element_shell_effective_plastic_strain] = {} - _ = array_requirements[ArrayType.element_shell_effective_plastic_strain] - # Radioss outputs needed to comptute Dyna output - _["dependents"] = ["element_shell_plastic_strain_upper", 'element_shell_plastic_strain_lower'] - _["shape"] = (1, n_shell, nip_shell) - _["convert"] = convert.element_shell_effective_plastic_strain - _["tracker"] = shell_ids_tracker - _["additional"] = [nip_shell] - - flag = "PARTS" + if n_solids > 0: + database_extent_binary[flag][0] += [ArrayType.element_solid_strain] + + # Dyna output + array_requirements[ArrayType.element_solid_strain] = {} + _ = array_requirements[ArrayType.element_solid_strain] + _["dependents"] = ["element_solid_strain"] + _["shape"] = (1, n_solids, nip_solid, 6) + _["convert"] = convert.element_solid_strain + _["tracker"] = solid_ids_tracker + _["additional"] = [nip_solid] + - database_extent_binary[flag] = {} - _ = database_extent_binary[flag] # Part masses for shells if rr.raw_header["nbEFunc"] > 0: + + flag = "PARTS" + + database_extent_binary[flag] = {} + _ = database_extent_binary[flag] database_extent_binary[flag][0] = [] database_extent_binary[flag][0] = _[0] + [ArrayType.part_mass] From 28878b5be375056ccefe23bacb4df51af37cc6ca Mon Sep 17 00:00:00 2001 From: Alexandre Hazo Date: Tue, 10 Jun 2025 10:05:18 +0200 Subject: [PATCH 3/3] fix rotate_tensor method --- vortex_radioss/animtod3plot/Anim_to_D3plot.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/vortex_radioss/animtod3plot/Anim_to_D3plot.py b/vortex_radioss/animtod3plot/Anim_to_D3plot.py index 76ed9d0..9af4c01 100644 --- a/vortex_radioss/animtod3plot/Anim_to_D3plot.py +++ b/vortex_radioss/animtod3plot/Anim_to_D3plot.py @@ -107,14 +107,13 @@ def rotate_tensor(vecs , R): Rt = np.transpose(R, axes=(0, 2, 1)) S_global = np.matmul(np.matmul(R, S_local), Rt) - # Extract the components and set the out-of-plane components to zero result = np.stack([ S_global[:, 0, 0], # ox S_global[:, 1, 1], # oy - np.zeros(n), # oz = 0 (out-of-plane) + S_global[:, 2, 2], # oz S_global[:, 0, 1], # oxy - np.zeros(n), # oyz = 0 (out-of-plane) - np.zeros(n) # oxz = 0 (out-of-plane) + S_global[:, 1, 2], # oyz + S_global[:, 0, 2], # oxz ], axis=-1) return result