utilities.py

import numpy as np
from scipy.interpolate import interp1d
from pylab import plt

def process_rpt(fname):
    with open(fname, 'rb') as f:
        lines = f.readlines()
    lines_cleaned = []
    for line in lines:
        # Strip leading and trailing whitespace
        line_cleaned = line.lstrip().rstrip()
        # Ignore blank lines
        if len(line_cleaned):
            # Check if first element is a float
            line_cleaned = line_cleaned.split()
            try:
                float(line_cleaned[0])
                lines_cleaned.append(line_cleaned)
            except:
                pass
    tmp = np.zeros((len(lines_cleaned),len(lines_cleaned[0])))
    for i,line in enumerate(lines_cleaned):
        tmp[i,:] = np.array(map(float, line))
    return tmp

class hot_spot_stress(object):
    def __init__(self, abaqus_xy_report_file, t, order=2):
        fname = abaqus_xy_report_file
        self.s_max_path = process_rpt(fname)
        self.t = t
        self.order = order
        self.x = self.s_max_path[:,0]
        self.s_max_range = self.s_max_path[:,1]

    def calculate(self):
        self.trialX = np.linspace(self.x[0], self.x[-1], 100)
        self.fitted_interp_2nd = np.polyfit(self.x, self.s_max_range, self.order)[::-1]
        self.second_order_y = np.zeros(len(self.trialX))
        for i in range(len(self.fitted_interp_2nd)):
            self.second_order_y += self.fitted_interp_2nd[i]*self.trialX**i
        f = interp1d(self.trialX, self.second_order_y)
        if self.x[-1] < 1.5*self.t:
            end = self.x[-1]
        else:
            end = 1.5*self.t
        interp_x = np.array([0.5*self.t,end])
        interp_y = f(interp_x)
        self.fitted_interp = np.polyfit(interp_x, interp_y, 1)[::-1]
        self.straight_y = np.zeros(len(self.trialX))
        for i in range(len(self.fitted_interp)):
            self.straight_y += self.fitted_interp[i]*self.trialX**i

    def plot(self):
        plt.figure()
        plt.plot(self.x,
                 self.s_max_range,
                 label='Data',
                 marker='o')
        plt.plot(self.trialX,
                 self.straight_y,
                 label = 'Linear extrapolation')
        plt.plot(self.trialX,
                 self.second_order_y,
                 label = 'Polynomial order {} fit'.format(self.order))
        plt.legend()
        plt.ymin=0.0

    def plot_raw(self):
        plt.figure()
        plt.plot(self.x,
                 self.s_max_range,
                 label='Data',
                 marker='o')
        plt.legend()
        plt.ymin=0.0

    def hot_spot(self, z):
        root_interp_stress = 0.0
        for i in range(len(self.fitted_interp)):
            root_interp_stress += self.fitted_interp[i]*z**i
        return root_interp_stress