Possible Mass Visualization

[Vachan-Vanian]

Hi @yexiang1992, I would like to know if the opstool could support mass visualization features via pyvista.

[yexiang92]

Hi @Vachan-Vanian , could you explain it in more detail?

[Vachan-Vanian]

Hi @yexiang1992,
I'm looking for a way to visually represent nodal mass distributions in the 3D model view, which would make it much easier to verify mass assignments and identify potential modeling issues. Specifically:

1. Mass representation as scaled points/spheres:
   ◦ Each node with mass would display a sphere with size proportional to its total translational mass (mx+my+mz) as in many cases the rotational is neglected. But it would be great if we could taking into account rotational too.
   ◦ This would provide an immediate visual cue about mass distribution throughout the structure
2. Color-coding scheme:
   ◦ Green spheres: Nodes where all 6 DOF values (mx, my, mz, rotx, roty, rotz) are identical
   ◦ Orange spheres: Nodes where the values are all zeros or are all none zero and not same values are identical
   ◦ Red spheres: Nodes with mixed assignment for same category (some zero, some non-zero values)
3. Mass value labels:
   ◦ Labels showing the actual mass values in two rows:
   ▪ Row 1: mx my mz (Translational)
   ▪ Row 2: rotx roty rotz (Rotational)
   ◦ These could be toggled on/off as needed to avoid visual clutter
4. Implementation using PyVista:
   ◦ Since opstool already utilizes PyVista for 3D visualization, these mass indicators could be integrated as additional actors in the same scene
   ◦ PyVista's ability to handle point glyphs with varying size and color would be perfect for this application
   This feature would be extremely helpful for model debugging and verification, especially for large structural models with complex mass distributions or when setting up dynamic analyses.

[yexiang92]

@Vachan-Vanian Hi! This is a great idea! However, there are some challenges to consider. The overall mass matrix in OpenSees is built from contributions at the nodes, elements, and within the materials themselves. Currently, we can only visualize mass applied directly at the nodes via nodeMass. If element or material properties contribute to the mass, this visualization will be incomplete. We need to explore how to account for these contributions to provide a comprehensive representation of the mass distribution.

[Vachan-Vanian]

@yexiang1992 hello! The most direct approach to overcome these challenges is to patch the OpenSees commands related to our goal. I provide below some draft code for a simple 2D problem where we manually create the matrix of masses based on the OpenSees commands that the user used, and then compare with the actual mass matrix from OpenSees to ensure that our calculation is correct.
If you agree with this approach, we can extend it to account for different model configurations and element types, materials, operations etc. also try to follow opstool api logic!
The logger could potentially used and for other parts of the codebase where it is not easy to retrieve information from OpenSees natively.

import inspect
from functools import wraps
import openseespy.opensees as ops
import numpy as np

# Dictionary to store mass-related data extracted from commands
mass_data = {
    'nodes': {},
    'elements': {},
    'nodal_masses': {},
    'model_info': {'ndm': 2, 'ndf': 3}  # Default values for the draft code
}

# List of commands to track
commands_to_include = ['node', 'element', 'mass', 'model']

def log_command(func):
    @wraps(func)
    def wrapper(*args, **kwargs):
        # Get the function name
        func_name = func.__name__
        
        # Only process if the function name is in the inclusion list
        if func_name in commands_to_include:

            if func_name == 'model':
                # Update model dimensions
                for i in range(len(args)):
                    if args[i] == '-ndm' and i+1 < len(args):
                        mass_data['model_info']['ndm'] = args[i+1]
                    elif args[i] == '-ndf' and i+1 < len(args):
                        mass_data['model_info']['ndf'] = args[i+1]

            elif func_name == 'node':
                # Node command: ops.node(tag, x, y, z, [options])
                node_tag = args[0]
                ndm = mass_data['model_info']['ndm']
                coords = args[1:1+ndm]  # Extract x,y,z based on ndm
                mass_data['nodes'][node_tag] = coords
                
                # Check for mass in node command
                for i in range(len(args)):
                    if args[i] == '-mass' and i+1 < len(args):
                        ndf = mass_data['model_info']['ndf']
                        # Get mass values based on the number specified and ndf
                        mass_values = []
                        for j in range(min(ndf, len(args) - (i+1))):
                            mass_values.append(args[i+1+j])
                        
                        # Pad with zeros if needed
                        while len(mass_values) < ndf:
                            mass_values.append(0.0)
                        
                        # Set nodal mass entry (replace any existing values)
                        mass_data['nodal_masses'][node_tag] = mass_values
            
            # Extract mass information from separate mass command
            elif func_name == 'mass':
                # Mass command: ops.mass(tag, mx, my, mz, ...) depending on ndf
                node_tag = args[0]
                ndf = mass_data['model_info']['ndf']
                mass_values = list(args[1:1+ndf])  # Extract masses based on ndf
                
                # Pad with zeros if needed
                while len(mass_values) < ndf:
                    mass_values.append(0.0)
                
                # Replace existing nodal mass entry
                mass_data['nodal_masses'][node_tag] = mass_values
            
            # Extract element information
            elif func_name == 'element':
                # General approach for any element type
                element_type = args[0]
                element_tag = args[1]
                
                # Create basic element info
                element_info = {
                    'type': element_type,
                    'nodes': [],
                    'mass_per_length': 0.0,
                    'is_consistent': False
                }
                
                # For most elements, the first two arguments after tag are nodes
                # But this can vary based on element type
                if len(args) > 3:
                    # Simple approach: assume first few args are nodes
                    # For more complex elements, this would need to be refined
                    num_nodes = 2  # Default assumption
                    element_info['nodes'] = list(args[2:2+num_nodes])
                
                # Search for mass-related keywords
                for i in range(2, len(args)):
                    if args[i] == '-mass' and i+1 < len(args):
                        element_info['mass_per_length'] = args[i+1]
                    elif args[i] == '-cMass':
                        element_info['is_consistent'] = True
                
                # Store the element info
                mass_data['elements'][element_tag] = element_info
        
        # Call the original function
        return func(*args, **kwargs)
    
    return wrapper

# Patch the OpenSees module
for name, obj in inspect.getmembers(ops):
    if inspect.isfunction(obj) or inspect.isbuiltin(obj):
        setattr(ops, name, log_command(obj))

print("OpenSees commands will be stored in memory for later calculations")

def calculate_manual_mass_matrix():
    """Calculate mass matrix manually based on logged commands"""
    # Calculate element lengths
    for tag, elem in mass_data['elements'].items():
        node1, node2 = elem['nodes']
        x1, y1 = mass_data['nodes'][node1]
        x2, y2 = mass_data['nodes'][node2]
        length = ((x2-x1)**2 + (y2-y1)**2)**0.5
        elem['length'] = length
    
    # Calculate mass matrix
    num_nodes = len(mass_data['nodes'])
    ndf = 3  # 3 DOFs per node in 2D
    total_size = num_nodes * ndf
    mass_matrix = np.zeros((total_size, total_size))
    
    # Add nodal masses
    for node_id, masses in mass_data['nodal_masses'].items():
        start_dof = (node_id-1) * ndf
        for i in range(len(masses)):
            mass_matrix[start_dof + i, start_dof + i] += masses[i]
    
    # Add element masses
    for elem_id, elem in mass_data['elements'].items():
        node1, node2 = elem['nodes']
        length = elem['length']
        mass_val = elem['mass_per_length']
        is_consistent = elem['is_consistent']
        
        # Total element mass
        elem_mass = mass_val * length
        
        # DOF indices
        i1 = (node1-1) * ndf
        i2 = (node2-1) * ndf
        
        if is_consistent:
            # Consistent mass formulation
            # Axial terms (u1, u2)
            mass_matrix[i1, i1] += elem_mass/3         # u1-u1
            mass_matrix[i1, i2] += elem_mass/6         # u1-u2
            mass_matrix[i2, i1] += elem_mass/6         # u2-u1
            mass_matrix[i2, i2] += elem_mass/3         # u2-u2
            
            # Transverse terms (v1, v2, θ1, θ2)
            mass_matrix[i1+1, i1+1] += elem_mass*13/35     # v1-v1
            mass_matrix[i1+1, i1+2] += elem_mass*11*length/210  # v1-θ1
            mass_matrix[i1+1, i2+1] += elem_mass*9/70      # v1-v2
            mass_matrix[i1+1, i2+2] += -elem_mass*13*length/420  # v1-θ2
            
            mass_matrix[i1+2, i1+1] += elem_mass*11*length/210  # θ1-v1
            mass_matrix[i1+2, i1+2] += elem_mass*length*length/105  # θ1-θ1
            mass_matrix[i1+2, i2+1] += elem_mass*13*length/420  # θ1-v2
            mass_matrix[i1+2, i2+2] += -elem_mass*length*length/140  # θ1-θ2
            
            mass_matrix[i2+1, i1+1] += elem_mass*9/70      # v2-v1
            mass_matrix[i2+1, i1+2] += elem_mass*13*length/420  # v2-θ1
            mass_matrix[i2+1, i2+1] += elem_mass*13/35     # v2-v2
            mass_matrix[i2+1, i2+2] += -elem_mass*11*length/210  # v2-θ2
            
            mass_matrix[i2+2, i1+1] += -elem_mass*13*length/420  # θ2-v1
            mass_matrix[i2+2, i1+2] += -elem_mass*length*length/140  # θ2-θ1
            mass_matrix[i2+2, i2+1] += -elem_mass*11*length/210  # θ2-v2
            mass_matrix[i2+2, i2+2] += elem_mass*length*length/105  # θ2-θ2
            
        else:
            # Lumped mass formulation
            mass_matrix[i1, i1] += elem_mass/2      # u1
            mass_matrix[i1+1, i1+1] += elem_mass/2  # v1
            mass_matrix[i2, i2] += elem_mass/2      # u2
            mass_matrix[i2+1, i2+1] += elem_mass/2  # v2
            # No mass for rotational DOFs
    
    return mass_matrix

def calculate_total_mass(matrix):
    """Calculate total mass using rigid body modes"""
    ndf = 3
    num_nodes = matrix.shape[0] // ndf
    
    # Create rigid body mode vectors
    x_mode = np.zeros(matrix.shape[0])
    y_mode = np.zeros(matrix.shape[0])
    
    for i in range(num_nodes):
        x_mode[i*ndf] = 1      # X-translation
        y_mode[i*ndf + 1] = 1  # Y-translation
    
    x_mass = x_mode @ matrix @ x_mode
    y_mass = y_mode @ matrix @ y_mode
    
    print(f"Total mass (X-direction): {x_mass:.2f} kg")
    print(f"Total mass (Y-direction): {y_mass:.2f} kg")
    
    return (x_mass + y_mass) / 2



# Simple exmaple
ops.wipe()
ops.model('Basic', '-ndm', 2, '-ndf', 3)

# Define nodes
ops.node(1, 0.0, 0.0, '-mass', 45, 70, 10) # add mass during node creation
ops.node(2, 5.0, 0.0)

# Define properties
A = 0.02
E = 2.1e11
Iz = 8.33e-6
mass_per_length = 157  # kg/m (ρA = 7850 * 0.02)

# Add nodal masses
# ops.mass(1, 25, 90, 0)  # mass in x, y directions and  rotational mass
ops.mass(2, 10, 75, 0.0)  # mass in x, y directions, no rotational mass

# Transformation
ops.geomTransf('Linear', 1)

# Define beam element with lumped mass
ops.element('elasticBeamColumn', 1, 1, 2, A, E, Iz, 1, '-mass', mass_per_length, '-cMass') # Consistent Mass
ops.element('elasticBeamColumn', 2, 2, 1, A, E, Iz, 1, '-mass', mass_per_length*2) # Lumed Mass

# Calculate manual mass matrix
manual_matrix = calculate_manual_mass_matrix()
print("\nManually calculated mass matrix from logged commands:")
print(np.array2string(manual_matrix, precision=2))

# Calculate total mass
total_mass = calculate_total_mass(manual_matrix)
print(f"Average total mass: {total_mass:.2f} kg")

# Compare with OpenSeesPy mass matrix
print("\nComparing with OpenSeesPy mass matrix:")
ops.wipeAnalysis()
ops.system('FullGeneral')
ops.numberer('Plain')
ops.constraints('Transformation')
ops.algorithm('Newton')
ops.integrator('Newmark', 0.5, 0.25)
ops.analysis('Transient')
ops.integrator('GimmeMCK', 1.0, 0.0, 0.0)  # Only mass
ops.analyze(1, 0.0)

# Get OpenSeesPy mass matrix
opensees_matrix = np.array(ops.printA('-ret')).reshape(6, 6)
print("OpenSeesPy mass matrix:")
print(np.array2string(opensees_matrix, precision=2))

# Calculate difference
max_diff = np.max(np.abs(manual_matrix - opensees_matrix))
print(f"Maximum difference between manual and OpenSeesPy: {max_diff:.6f}")

---

OUTPUT CONSOLE

OpenSees commands will be stored in memory for later calculations

Manually calculated mass matrix from logged commands:
[[1091.67    0.      0.    130.83    0.      0.  ]
 [   0.   1146.57  205.6     0.    100.93 -121.49]
 [   0.    205.6   196.9     0.    121.49 -140.18]
 [ 130.83    0.      0.   1056.67    0.      0.  ]
 [   0.    100.93  121.49    0.   1151.57 -205.6 ]
 [   0.   -121.49 -140.18    0.   -205.6   186.9 ]]
Total mass (X-direction): 2410.00 kg
Total mass (Y-direction): 2500.00 kg
Average total mass: 2455.00 kg

Comparing with OpenSeesPy mass matrix:
OpenSeesPy mass matrix:
[[1091.67    0.      0.    130.83    0.      0.  ]
 [   0.   1146.57  205.6     0.    100.93 -121.49]
 [   0.    205.6   196.9     0.    121.49 -140.18]
 [ 130.83    0.      0.   1056.67    0.      0.  ]
 [   0.    100.93  121.49    0.   1151.57 -205.6 ]
 [   0.   -121.49 -140.18    0.   -205.6   186.9 ]]
Maximum difference between manual and OpenSeesPy: 0.000000

looking forward for your comments!

P.S. I have reviewed the source code of the opstool, and I can tell you have already implemented the node part very well.
For the element mass, although I think users need to specifically define them with opstool commands and not with general OpenSees commands?

[yexiang92]

@Vachan-Vanian Hi there,

This looks like a really good idea—let me take some time to think it over. Besides mass, there might be more great things we can do. Thank you for your contribution; it seems like a fairly general approach:

Extracting all relevant data entered by users through openseespy commands, including (but not limited to) mass, element properties, and section properties.

Best regards!

[Vachan-Vanian]

Hi @yexiang1992,
Thank you for considering my suggestion! I'm glad you see the potential in this approach for extracting and visualizing mass data. Indeed this method could be extended beyond just mass visualization to other properties as well.