Skip to content

Implement base transformation from local to global for stress and strain tensors (shell elements) #19

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
alexandrehazo wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Implement base transformation from local to global for stress and strain tensors (shell elements) #19

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
478145d
Implement base transformation from local to global for stress and str…
May 5, 2025
44368a8
Vectorization of local axes computation for shell elements
May 9, 2025
28878b5
fix rotate_tensor method
Jun 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
269 changes: 232 additions & 37 deletions vortex_radioss/animtod3plot/Anim_to_D3plot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -75,6 +75,48 @@ def sum2Scalars(*data):
print(data1)
print(data2)
return sum(data1) + sum(data2)

@staticmethod
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):
"""
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)

result = np.stack([
S_global[:, 0, 0], # ox
S_global[:, 1, 1], # oy
S_global[:, 2, 2], # oz
S_global[:, 0, 1], # oxy
S_global[:, 1, 2], # oyz
S_global[:, 0, 2], # oxz
], axis=-1)

return result

@staticmethod
def element_shell_internal_energy(*data):
Expand All @@ -89,29 +131,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
Expand All @@ -132,8 +191,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__(
Expand Down Expand Up @@ -258,9 +336,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:

Expand Down Expand Up @@ -546,22 +631,97 @@ 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]

# [0] are essential arrays and need creating even if no data available

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_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]



# 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"]

# 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)

# 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] = {}
Expand Down Expand Up @@ -601,8 +761,8 @@ 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_effective_plastic_strain] = {}
_ = array_requirements[ArrayType.element_shell_effective_plastic_strain]
Expand All @@ -612,14 +772,49 @@ def A_2_D(self,file_stem, use_shell_mask, use_solid_mask, use_beam_mask, silent,
_["convert"] = convert.element_shell_effective_plastic_strain
_["tracker"] = shell_ids_tracker
_["additional"] = [nip_shell]

flag = "PARTS"



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]
# 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]

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]



# 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]
Expand Down