cpress export for Abaqus

Abaqus Script.py

from abaqus import *
from abaqusConstants import *
from odbAccess import *

def extract_node_data(odb_path, step_name, frame_number=-1):
    """
    Extract node coordinates and CPRESS field from an Abaqus ODB file
    
    Parameters:
    -----------
    odb_path : str
        Path to the ODB file
    step_name : str
        Name of the step to extract data from
    frame_number : int
        Frame number to extract (default: -1, last frame)
    
    Returns:
    --------
    node_coords : list
        List of node coordinates [[x,y,z], ...]
    cpress : list
        List of CPRESS values for each node (None for nodes without CPRESS)
    nodes_with_cpress : list
        List of node labels that have CPRESS values
    """
    
    # Open the ODB file
    odb = openOdb(path=odb_path)
    
    try:
        # Get the specified step
        step = odb.steps[step_name]
        
        # Get the last frame if frame_number is -1
        if frame_number == -1:
            frame = step.frames[-1]
        else:
            frame = step.frames[frame_number]
        
        # Get the node coordinates
        assembly = odb.rootAssembly
        instance = assembly.instances['PART-1-1']  # Modify instance name if needed
        nodes = instance.nodes
        
        # Initialize lists for storing data
        node_coords = []
        node_labels = []  # To store original node labels
        
        for node in nodes:
            coords = [node.coordinates[0], 
                     node.coordinates[1], 
                     node.coordinates[2]]
            node_coords.append(coords)
            node_labels.append(node.label)
        
        # Get CPRESS field
        field_output = frame.fieldOutputs['CPRESS']
        num_nodes = len(nodes)
        cpress_values = [None] * num_nodes  # Initialize with None
        nodes_with_cpress = []  # To store nodes that actually have CPRESS values
        
        # Create a mapping from node label to index
        label_to_index = {}
        for i, label in enumerate(node_labels):
            label_to_index[label] = i
        
        # Process CPRESS values
        for value in field_output.values:
            node_label = value.nodeLabel
            if node_label in label_to_index:
                idx = label_to_index[node_label]
                cpress_values[idx] = value.data
                nodes_with_cpress.append(node_label)
        
        return node_coords, cpress_values, nodes_with_cpress
        
    finally:
        # Close the ODB
        odb.close()

def save_to_file(node_coords, cpress, nodes_with_cpress, output_file='output_data.txt'):
    """
    Save the extracted data to text files
    """
    # Save all nodes with their coordinates and CPRESS (if available)
    with open(output_file, 'w') as f:
        f.write("Node_Index Node_X Node_Y Node_Z CPRESS\n")
        for i in xrange(len(node_coords)):
            cpress_val = cpress[i]
            if cpress_val is None:
                cpress_str = "None"
            else:
                cpress_str = "%.6f" % cpress_val
            
            f.write("%d %f %f %f %s\n" % (
                i+1,
                node_coords[i][0],
                node_coords[i][1],
                node_coords[i][2],
                cpress_str
            ))
    
    # Save a separate file with only nodes that have CPRESS values
    contact_file = output_file.replace('.txt', '_contact_only.txt')
    with open(contact_file, 'w') as f:
        f.write("Node_Label Node_X Node_Y Node_Z CPRESS\n")
        for node_label in nodes_with_cpress:
            idx = node_label - 1  # Assuming node labels start from 1
            f.write("%d %f %f %f %f\n" % (
                node_label,
                node_coords[idx][0],
                node_coords[idx][1],
                node_coords[idx][2],
                cpress[idx]
            ))

def main():
    """
    Main function to run the extraction
    """
    # Configuration
    odb_path = 'Cpress_02.odb'  # Change this
    step_name = 'Step-1'  # Change this to match your step name
    output_file = 'output_data.txt'
    
    try:
        print('Starting data extraction from ODB...')
        node_coords, cpress, nodes_with_cpress = extract_node_data(odb_path, step_name)
        
        print('Saving data to files...')
        save_to_file(node_coords, cpress, nodes_with_cpress, output_file)
        
        print('Successfully extracted data to:')
        print('- All nodes: %s' % output_file)
        print('- Contact nodes only: %s' % output_file.replace('.txt', '_contact_only.txt'))
        print('Total number of nodes: %d' % len(node_coords))
        print('Number of nodes with CPRESS: %d' % len(nodes_with_cpress))
        
    except Exception, e:  # Python 2 exception syntax
        print('Error occurred: %s' % str(e))

# Run the script
if __name__ == "__main__":
    main()

Visualization.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from scipy.interpolate import griddata, LinearNDInterpolator
from pathlib import Path


def read_contact_data(filepath):
    """
    Read contact data from the contact-only output file
    """
    data = np.loadtxt(filepath, skiprows=1)
    return data[:, 1:5]  # x, y, z, cpress


def cylindrical_to_unwrapped(x, y, z, cpress, axis='z'):
    """
    Convert cylindrical coordinates to unwrapped 2D coordinates
    """
    if axis == 'z':
        r = np.sqrt(x ** 2 + y ** 2)
        theta = np.arctan2(y, x)
        height = z
    elif axis == 'y':
        r = np.sqrt(x ** 2 + z ** 2)
        theta = np.arctan2(z, x)
        height = y
    else:  # x axis
        r = np.sqrt(y ** 2 + z ** 2)
        theta = np.arctan2(z, y)
        height = x

    # Convert theta from [-pi, pi] to [0, 2pi]
    theta = np.where(theta < 0, theta + 2 * np.pi, theta)

    return theta, height, r, cpress


def handle_periodicity(theta_deg, height, cpress, width=10):
    """
    Handle periodicity by adding duplicate points at boundaries
    """
    # Add points at 360° that duplicate those at 0°
    mask_left = theta_deg < width
    mask_right = theta_deg > (360 - width)

    # Add points beyond both boundaries
    theta_extended = np.concatenate([
        theta_deg[mask_right] - 360,  # Points before 0°
        theta_deg,
        theta_deg[mask_left] + 360  # Points after 360°
    ])

    height_extended = np.concatenate([
        height[mask_right],
        height,
        height[mask_left]
    ])

    cpress_extended = np.concatenate([
        cpress[mask_right],
        cpress,
        cpress[mask_left]
    ])

    return theta_extended, height_extended, cpress_extended


def plot_unwrapped_contact(data, axis='z', num_grid_points=200,
                           cmap='viridis', title=None, show_colorbar=True):
    """
    Create an unwrapped 2D plot of the cylindrical contact surface
    """
    # Extract coordinates and pressure
    x, y, z, cpress = data.T

    # Convert to unwrapped coordinates
    theta, height, radius, cpress = cylindrical_to_unwrapped(x, y, z, cpress, axis)

    # Convert theta to degrees for plotting
    theta_deg = np.degrees(theta)

    # Handle periodicity by adding points across the boundary
    theta_ext, height_ext, cpress_ext = handle_periodicity(theta_deg, height, cpress)

    # Add small amount of noise to avoid exact duplicates
    theta_ext += np.random.normal(0, 0.01, theta_ext.shape)
    height_ext += np.random.normal(0, 0.01, height_ext.shape)

    # Create regular grid for interpolation
    theta_grid = np.linspace(0, 360, num_grid_points)
    height_grid = np.linspace(height.min(), height.max(), num_grid_points)
    THETA, HEIGHT = np.meshgrid(theta_grid, height_grid)

    # Use LinearNDInterpolator instead of Rbf
    interp = LinearNDInterpolator(list(zip(theta_ext, height_ext)), cpress_ext, fill_value=np.nan)
    pressure_grid = interp(THETA, HEIGHT)

    # Create mask for valid pressure points
    mask = np.isfinite(pressure_grid)

    # Apply mask and clip values to original range
    pressure_grid = np.where(mask, pressure_grid, np.nan)
    pressure_grid = np.clip(pressure_grid, np.min(cpress), np.max(cpress))

    # Create figure
    plt.figure(figsize=(12, 8))

    # Plot the interpolated pressure
    im = plt.pcolormesh(THETA, HEIGHT, pressure_grid,
                        shading='gouraud',
                        cmap=cmap,
                        vmin=np.min(cpress),
                        vmax=np.max(cpress))

    if show_colorbar:
        cbar = plt.colorbar(im)
        cbar.set_label('Contact Pressure', rotation=270, labelpad=15)

    plt.xlabel('Angular Position (degrees)')
    plt.ylabel(f'Position along {axis}-axis')

    if title:
        plt.title(title)

    # Set the x-axis ticks to show degrees
    plt.xticks(np.linspace(0, 360, 9))

    plt.grid(True, linestyle='--', alpha=0.3)

    # Add some statistical information
    stats_text = f'Pressure Range: {np.min(cpress):.2f} to {np.max(cpress):.2f}\n'
    stats_text += f'Number of contact points: {len(cpress)}'
    plt.text(0.02, 0.98, stats_text,
             transform=plt.gca().transAxes,
             verticalalignment='top',
             bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

    return plt.gcf()


def main():
    # Configuration
    input_file = 'output_data_contact_only.txt'
    axis = 'z'

    try:
        # Read the data
        print("Reading contact data...")
        data = read_contact_data(input_file)

        print(f"Contact pressure range: {np.min(data[:, 3]):.2f} to {np.max(data[:, 3]):.2f}")

        # Create the plot
        print("Creating visualization...")
        fig = plot_unwrapped_contact(
            data,
            axis=axis,
            title=f'Unwrapped Cylindrical Contact Surface (Axis: {axis})',
            num_grid_points=200,  # Increased for smoother visualization
            cmap='viridis'
        )

        # Save the plot
        output_file = 'contact_visualization.png'
        fig.savefig(output_file, dpi=300, bbox_inches='tight')
        print(f"Visualization saved to: {output_file}")

        # Show the plot
        plt.show()

    except Exception as e:
        print(f"An error occurred: {str(e)}")


if __name__ == "__main__":
    main()