Make connector behavior and generate plots of stiffness

MkConnectors.py

import os
import sys
import odbAccess
from glob import glob
import re
from collections import defaultdict
import numpy as NP
import subprocess

def open_odb(odbPath):
    base, ext = os.path.splitext(odbPath)
    odbPath = base + '.odb'
    new_odbPath = None
    if odbAccess.isUpgradeRequiredForOdb(upgradeRequiredOdbPath=odbPath):
        print('odb %s needs upgrading' % (odbPath,))
        path,file_name = os.path.split(odbPath)
        file_name = base + "_upgraded.odb"
        new_odbPath = os.path.join(path,file_name)
        odbAccess.upgradeOdb(existingOdbPath=odbPath, upgradedOdbPath=new_odbPath)
        odbPath = new_odbPath
    odb = odbAccess.openOdb(path=odbPath, readOnly=True)
    return odb

def flatten_stupidly(l):
    out = []
    for item in l:
        try:
            len(item)
            out.extend(flatten_stupidly(item))
        except:
            out.append(item)
    return out


class VarRegEx(object):
    def __init__(self, regex, x_key, y_key):
        self.regex = regex
        self.x_key = x_key
        self.y_key = y_key

class NodeHistoryData(object):
    def __init__(self, odb, nset_name):
        self.odb = odb
        self.nset_name = nset_name
        self.history_rgn_keys = self.mk_nodal_history_rgn_keys()


    def extract(self, x_key, y_key):
        for history_rgn_key in self.history_rgn_keys:
            times, x_vals = self.get_history_data(history_rgn_key, x_key)
            times, y_vals = self.get_history_data(history_rgn_key, y_key)
            self.data = [(x,y) for x,y in zip(x_vals,y_vals)]

    def mk_nodal_history_rgn_keys(self):
        '''
        Need to make a key like "Node WELDFLAW-1.109"
        '''
        keys = []
        self.nodes = flatten_stupidly(self.odb.rootAssembly.nodeSets[self.nset_name].nodes)
        for node in self.nodes:
            if node.instanceName == '__unknown__':
                keys.append("Node ASSEMBLY." + str(node.label))
            else:
                keys.append("Node " + node.instanceName + "." + str(node.label))
        return keys

    def get_history_data(self, hist_rgn_key, result_key):
        times = []
        vals = []
        for step_name in odb.steps.keys():
            step = odb.steps[step_name]
            hist_rgn = step.historyRegions[hist_rgn_key]
            data = hist_rgn.historyOutputs[result_key].data
            times.extend([t + step.totalTime for t, val in data])
            vals.extend([val for t, val in data])
        return times, vals

def reflect(x_y_pairs,scale=1.0):
    x_y_array = NP.array([[x,scale*y] for x,y in x_y_pairs])
    x_y_array = NP.append(x_y_array,NP.array([[-x,-scale*y] for x,y in x_y_pairs[1:]]),axis=0)
    # Sort by x (first column)
    # See http://stackoverflow.com/questions/2828059/sorting-arrays-in-numpy-by-column
    x_y_array = x_y_array[x_y_array[:,0].argsort(),:]
    return x_y_array[:,0], x_y_array[:,1]

def combine(x_y_pairs_1, x_y_pairs_2, scale=1.0):
    x_y_array = NP.array([[x,scale*y] for x,y in x_y_pairs_1])
    x_y_array = NP.append(x_y_array,NP.array([[x,scale*y] for x,y in x_y_pairs_2[1:]]),axis=0)
    # Sort by x (first column)
    # See http://stackoverflow.com/questions/2828059/sorting-arrays-in-numpy-by-column
    x_y_array = x_y_array[x_y_array[:,0].argsort(),:]
    return x_y_array[:,0], x_y_array[:,1]

def mk_conn_behavior(behavior_name, node_histories):
    f = open(behavior_name + '.inp', 'wb')
    f.write('*Connector Behavior, name=%s, extrapolation=LINEAR\n' % (behavior_name,))
    # Component 1 is shear
    f.write('*Connector Elasticity, nonlinear, component=1\n')
    x,y = reflect(node_histories['Shear'][0].data, scale=2.0)
    for _x,_y in zip(x,y):
        f.write('%s,%s\n' % (_y,_x))
    xy_secant = NP.column_stack((x,y/x))
    xy_secant = xy_secant[~NP.isnan(xy_secant).any(1)]

    make_asy_plot(behavior_name+'Shear.asy',
                  x[:-1],
                  NP.diff(y)/NP.diff(x),
                  r"\mbox{Displacement} (in)",
                  r"\mbox{Tangent Stiffness} (lbf \cdot in)",
                  x2=xy_secant[:,0],
                  y2=xy_secant[:,1],
                  yaxis2=r"\mbox{Secant Stiffness} (lbf \cdot in)")
    # Component 2 is shear and = to 1
    f.write('*Connector Elasticity, nonlinear, component=2\n')
    for _x,_y in zip(x,y):
        f.write('%s,%s\n' % (_y,_x))
    # Component 3 is compression and tension
    f.write('*Connector Elasticity, nonlinear, component=3\n')
    x,y = combine(node_histories['Tension'][0].data, node_histories['Compression'][0].data,scale=2.0)
    xy_secant = NP.column_stack((x,y/x))
    xy_secant = xy_secant[~NP.isnan(xy_secant).any(1)]

    make_asy_plot(behavior_name+'Axial.asy',
                  x[:-1],
                  NP.diff(y)/NP.diff(x),
                  r"\mbox{Displacement} (in)",
                  r"\mbox{Tangent Stiffness} (lbf \cdot in)",
                  x2=xy_secant[:,0],
                  y2=xy_secant[:,1],
                  yaxis2=r"\mbox{Secant Stiffness} (lbf \cdot in)")
    for _x,_y in zip(x,y):
        f.write('%s,%s\n' % (_y,_x))
    # Component 4 is rotation
    f.write('*Connector Elasticity, nonlinear, component=4\n')
    x,y = reflect(node_histories['Rotation'][0].data, scale=2.0)
    for _x,_y in zip(x,y):
        f.write('%s,%s\n' % (_y,_x))
    xy_secant = NP.column_stack((x,y/x))
    xy_secant = xy_secant[~NP.isnan(xy_secant).any(1)]

    make_asy_plot(behavior_name+'Rotation.asy',
                  x[:-1],
                  NP.diff(y)/NP.diff(x),
                  r"\mbox{Rotation} (rad)",
                  r"\mbox{Tangent Stiffness} (\frac{lbf\cdot in}{rad})",
                  x2=xy_secant[:,0],
                  y2=xy_secant[:,1],
                  yaxis2=r"\mbox{Secant Stiffness} (\frac{lbf\cdot in}{rad})")

    # Component 5 is rotation = component 4
    f.write('*Connector Elasticity, nonlinear, component=5\n')
    for _x,_y in zip(x,y):
        f.write('%s,%s\n' % (_y,_x))
    # Component 6 is torsion. Scale = 1 since full model is used
    f.write('*Connector Elasticity, nonlinear, component=6\n')
    x,y = reflect(node_histories['Torsion'][0].data, scale=1.0)
    for _x,_y in zip(x,y):
        f.write('%s,%s\n' % (_y,_x))
    xy_secant = NP.column_stack((x,y/x))
    xy_secant = xy_secant[~NP.isnan(xy_secant).any(1)]

    make_asy_plot(behavior_name+'Torsion.asy',
                  x[:-1],
                  NP.diff(y)/NP.diff(x),
                  r"\mbox{Rotation} (rad)",
                  r"\mbox{Tangent Stiffness} (\frac{lbf\cdot in}{rad})",
                  x2=xy_secant[:,0],
                  y2=xy_secant[:,1],
                  yaxis2=r"\mbox{Secant Stiffness} (\frac{lbf\cdot in}{rad})")
    f.close()

    # Cocking stiffness
    x,y = reflect(node_histories['Cocking'][0].data, scale=2.0)
    xy_secant = NP.column_stack((x,y/x))
    xy_secant = xy_secant[~NP.isnan(xy_secant).any(1)]

    make_asy_plot(behavior_name+'Cocking.asy',
                  x[:-1],
                  NP.diff(y)/NP.diff(x),
                  r"\mbox{Rotation} (rad)",
                  r"\mbox{Tangent Stiffness} (\frac{lbf\cdot in}{rad})",
                  x2=xy_secant[:,0],
                  y2=xy_secant[:,1],
                  yaxis2=r"\mbox{Secant Stiffness} (\frac{lbf\cdot in}{rad})")

def make_asy_plot(name, x, y, xaxis, yaxis, x2=None, y2=None, yaxis2=None ):
    '''
    --col_x Y
    --col_y Y
    --caption_x "\mbox{Rotation} (rad)"
    --caption_y "\mbox{Secant Stiffness} (\frac{lbf\cdot in}{rad})"
     --col_x Y
     --col_y Y
     --caption_x "\mbox{Displacement} (in)"
     --caption_y "\mbox{Secant Stiffness} (lbf \cdot in)"
    make_asy_plot(args.output_asy, x[1:], y[1:], args.caption_x, args.caption_y)
    '''
    f = open(name, 'wb')
    f.write('import graph;\n')
    f.write('pen dotted=linetype(new real[] {0,4});\n')
    f.write('size(400,300,IgnoreAspect);\n')
    f.write('real[] x={%s};\n' % (','.join([str(val) for val in x]),))
    f.write('real[] y={%s};\n' % (','.join([str(val) for val in y]),))

    f.write('scale(xautoscale=true, yautoscale=true);\n')
    f.write('draw(graph(x,y),black);\n')
    f.write('xaxis("$%s$",BottomTop,LeftTicks(N=8, extend=true, n=5, ptick=dotted));\n' % (xaxis,))
    f.write('yaxis("$%s$",LeftRight,RightTicks(N=10,extend=true, n=5, ptick=dotted));\n' % (yaxis,))
    if x2 is not None and y2 is not None:
        f.write('real[] x2={%s};\n' % (','.join([str(val) for val in x2]),))
        f.write('real[] y2={%s};\n' % (','.join([str(val) for val in y2]),))
        f.write('picture secondary=secondaryY(new void(picture pic) {\n')
        f.write('scale(pic,Linear(true),Linear(true));\n')
        f.write('draw(pic,graph(pic,x2,y2),red);\n')
        f.write('yaxis(pic,"$%s$",Right,red,LeftTicks(N=10,begin=false,end=false));\n' % (yaxis2,))
        f.write('});\n')
        f.write('add(secondary);\n')

    f.close()
    proc = subprocess.Popen('asy -f pdf %s' % (name,),
                            shell=True,
                            stdout=subprocess.PIPE,
                            stderr=subprocess.STDOUT,
                            )

reg_exes = {'Cocking':VarRegEx(re.compile(r'_Cocking'), 'UR3', 'RM3'),
            'Compression':VarRegEx(re.compile(r'_Compression'), 'U2', 'RF2'),
            'Rotation':VarRegEx(re.compile(r'_Rotation'), 'UR3', 'RM3'),
            'Tension':VarRegEx(re.compile(r'_Tension'), 'U2', 'RF2'),
            'Torsion':VarRegEx(re.compile(r'_Torsion'), 'UR2', 'RM2'),
            'Shear':VarRegEx(re.compile(r'_X'), 'U1', 'RF1')}

def extract_var_key_from_odb(odb_name):
    for name, var_reg_ex in reg_exes.items():
        if var_reg_ex.regex.search(odb_name):
            return name, var_reg_ex.x_key, var_reg_ex.y_key


if __name__ == '__main__':
    nsets = ['TOP', ]

    behavior_name = sys.argv[1]
    odb_pattern = sys.argv[2]
    odbs = glob(odb_pattern+'*'+'.odb')
    node_histories = defaultdict(list)
    for odb in odbs:
        print('Working on %s' % (odb,))
        name, x_key, y_key = extract_var_key_from_odb(odb)
        base, ext = os.path.splitext(odb)
        odb = open_odb(odb)

        for nset_name in nsets:
            nh = NodeHistoryData(odb, nset_name)
            nh.extract(x_key, y_key)
            node_histories[name].append(nh)
        odb.close()
    mk_conn_behavior(behavior_name, node_histories)