'''
Created on 2013-04-24

@author: Scott Sibole
'''
import sys, string, gts, subprocess, os
import numpy as np

def main(cfg):
    ccontent = parse_cfg(cfg)
    blocks = []
    cnt = 0
    for c in ccontent:
        blocks.append(make_block(c,cnt))
        cnt += 1
    fid = open(string.replace(cfg,'.cfg','.tg'),'w')
    for b in blocks:
        for line in b:
            fid.write(line)
    fid.close()
            
    
def parse_cfg(cfg):
    fid = open(cfg,'r')
    lines = map(string.rstrip,fid.readlines())
    fid.close()
    ccontent = []
    cnt = 1
    for line in lines:
        dmy = line.split(',')
        try:
            ccontent.append([dmy[0],dmy[1],float(dmy[2]),cnt])
            cnt += 2
        except:
            print 'Warning: Could not interpret configuration line containing string:', dmy
    return ccontent

def make_block(c,blkcnt):
    nf = [] #list of gts filenames
    s = []  #list of surfaces
    bounds = [] #list of bounding boxes
    nelms = []  # number of elements to use [i,j,k]
    for i in xrange(2):
        #convert from stl to gts format
        nf.append(string.replace(c[i],'.stl','.gts'))
        subprocess.call('stl2gts < '+c[i]+' > '+nf[i], shell=True)
        #read in surface in gts format
        fid = open(nf[i],'r')
        s.append(gts.read(fid))
        fid.close()
        os.remove(nf[i])
        #get bounding boxes
        verts = s[i].vertices()
        a = np.zeros((len(verts),3),float)
        for j in xrange(len(verts)):
            a[j,0] = verts[j].x
            a[j,1] = verts[j].y
            a[j,2] = verts[j].z
        # calculate the mean
        mu = np.mean(a,axis=0)
        cv = np.dot((a-mu).T,(a-mu))/float(len(verts)-1)
        w, v = np.linalg.eig(cv)
    #               8----------7
    #              /          /|
    #             5---------6  |
    #             |         |  |
    #             |   4     | 3|
    #             |         | /
    #             |1_______2|/

        v = abs(v)
        if np.argmax(v[:,0]) == 0:
            vi = v[:,0]
            ri = w[0]
            if np.argmax(v[:,1]) == 1:
                vj = v[:,1]
                rj = w[1]
                vk = v[:,2]
                rk = w[2]
            else:
                vj = v[:,2]
                rj = w[2]
                vk = v[:,1]
                rk = w[1]
        elif np.argmax(v[:,0]) == 1:
            vj = v[:,0]
            rj = w[0]
            if np.argmax(v[:,1]) == 0:
                vi = v[:,1]
                ri = w[1]
                vk = v[:,2]
                rk = w[2]
            else:
                vi = v[:,2]
                ri = w[2]
                vk = v[:,1]
                rk = w[1]
        else:
            vk = v[:,0]
            rk = w[0]
            if np.argmax(v[:,1]) == 0:
                vi = v[:,1]
                ri = w[1]
                vj = v[:,2]
                rj = w[2]
            else:
                vi = v[:,2]
                ri = w[2]
                vj = v[:,1]
                rj = w[1]
        ri = np.sqrt(ri)
        rj = np.sqrt(rj)
        rk = np.sqrt(rk)
        c1 = mu-ri*vi-rj*vj-rk*vk
        c2 = mu+ri*vi-rj*vj-rk*vk
        c3 = mu+ri*vi+rj*vj-rk*vk
        c4 = mu-ri*vi+rj*vj-rk*vk
        c5 = mu-ri*vi-rj*vj+rk*vk
        c6 = mu+ri*vi-rj*vj+rk*vk
        c7 = mu+ri*vi+rj*vj+rk*vk
        c8 = mu-ri*vi+rj*vk+rk*vk
        bounds.append(np.vstack((c1,c2,c3,c4,c5,c6,c7,c8)))
        nelms.append([np.ceil(2*ri/c[2]),np.ceil(2*rj/c[2]),np.ceil(2*rk/c[2])])

    ep = []
    ec = []
    for i in xrange(3):
        ep.append(nelms[1][i]-nelms[0][i])
        ec.append(nelms[0][i])

    tg_lines = []
    tg_lines.append('\nc Cell and PCM '+str((c[3]+1)/2)+'\n')
    tg_lines.append('c Import Surfaces\n')
    tg_lines.append('sd '+str(blkcnt*2+1)+' stl '+c[0]+';\n')
    tg_lines.append('sd '+str(blkcnt*2+2)+' stl '+c[1]+';\n\nc Create Block\n')
    blk_str = 'block '
    for i in xrange(3):
        en = [1,int(1+ep[i]),int(1+2*ep[i]),int(1+ec[i]+2*ep[i]),int(1+ec[i]+3*ep[i]),int(1+ec[i]+4*ep[i])]
        blk_str += str(en[0])+' '+str(en[1])+' '+str(en[2])+' '+str(en[3])+' '+str(en[4])+' '+str(en[5])+'; '
    tg_lines.append(blk_str+'\n')
    for i in xrange(3):
        tg_lines.append('0 0 0 0 0 0;\n')
        
    tg_lines.append('dei 1 3 0 4 6;; 1 3 0 4 6;\n')
    tg_lines.append('dei 3 4; 1 3 0 4 6; 1 3 0 4 6;\n')
    tg_lines.append('dei 1 3 0 4 6; 1 3 0 4 6; 3 4;\n')
    tg_lines.append('\nc Reposition vertices\n')
    #pcm faces
    ijk1 = [[1,3,3],[6,3,3],[6,4,3],[1,4,3],[1,3,4],[6,3,4],[6,4,4],[1,4,4]]
    ijk2 = [[2,3,3],[5,3,3],[5,4,3],[2,4,3],[2,3,4],[5,3,4],[5,4,4],[2,4,4]]
    for i in xrange(8):
        tg_lines.append('pb '+string.join(map(str,ijk1[i]),' ')+' '+string.join(map(str,ijk1[i]),' ')+' xyz '+string.join(map(str,bounds[1][i,:]),' ')+'\n')
        tg_lines.append('pb '+string.join(map(str,ijk2[i]),' ')+' '+string.join(map(str,ijk2[i]),' ')+' xyz '+string.join(map(str,bounds[0][i,:]),' ')+'\n')
    ijk1 = [[3,1,3],[4,1,3],[4,6,3],[3,6,3],[3,1,4],[4,1,4],[4,6,4],[3,6,4]]
    ijk2 = [[3,2,3],[4,2,3],[4,5,3],[3,5,3],[3,2,4],[4,2,4],[4,5,4],[3,5,4]]
    for i in xrange(8):
        tg_lines.append('pb '+string.join(map(str,ijk1[i]),' ')+' '+string.join(map(str,ijk1[i]),' ')+' xyz '+string.join(map(str,bounds[1][i,:]),' ')+'\n')
        tg_lines.append('pb '+string.join(map(str,ijk2[i]),' ')+' '+string.join(map(str,ijk2[i]),' ')+' xyz '+string.join(map(str,bounds[0][i,:]),' ')+'\n')
    ijk1 = [[3,3,1],[4,3,1],[4,4,1],[3,4,1],[3,3,6],[4,3,6],[4,4,6],[3,4,6]]
    ijk2 = [[3,3,2],[4,3,2],[4,4,2],[3,4,2],[3,3,5],[4,3,5],[4,4,5],[3,4,5]]
    for i in xrange(8):
        tg_lines.append('pb '+string.join(map(str,ijk1[i]),' ')+' '+string.join(map(str,ijk1[i]),' ')+' xyz '+string.join(map(str,bounds[1][i,:]),' ')+'\n')
        tg_lines.append('pb '+string.join(map(str,ijk2[i]),' ')+' '+string.join(map(str,ijk2[i]),' ')+' xyz '+string.join(map(str,bounds[0][i,:]),' ')+'\n')
    #cell faces
    ijk1 = [[3,3,3],[4,3,3],[4,4,3],[3,4,3],[3,3,4],[4,3,4],[4,4,4],[3,4,4]]
    for i in xrange(8):
        tg_lines.append('pb '+string.join(map(str,ijk1[i]),' ')+' '+string.join(map(str,ijk1[i]),' ')+' xyz '+string.join(map(str,bounds[0][i,:]),' ')+'\n')
        
    '''
    i1 = 0
    i2 = 3
    for i in [1,6]:
        j1 = 1
        j2 = 1
        for j in [3,4]:
            k1 = 2
            k2 = 1
            for k in [3,4]:
                tg_lines.append('pb '+str(i)+' '+str(j)+' '+str(k)+' '+str(i)+' '+str(j)+' '+str(k)+' xyz '+str(bounds[1][i1,0])+' '+str(bounds[1][j1])+' '+str(bounds[1][k1])+'\n')
                tg_lines.append('pb '+str(i2)+' '+str(j2)+' '+str(k2)+' '+str(i2)+' '+str(j2)+' '+str(k2)+' xyz '+str(bounds[1][i1])+' '+str(bounds[1][j1])+' '+str(bounds[1][k1])+'\n')
                k1 = 5
                k2 = 6
            j1 = 4
            j2 = 6
        i1 = 3
        i2 = 4
        
    
    #cell faces
    i1 = 0
    i2 = 3
    for i in [2,5]:
        j1 = 1
        j2 = 2
        for j in [3,4]:
            k1 = 2
            k2 = 2
            for k in [3,4]:
                tg_lines.append('pb '+str(i)+' '+str(j)+' '+str(k)+' '+str(i)+' '+str(j)+' '+str(k)+' xyz '+str(bounds[1][i1])+' '+str(bounds[1][j1])+' '+str(bounds[1][k1])+'\n')
                tg_lines.append('pb '+str(i2)+' '+str(j2)+' '+str(k2)+' '+str(i2)+' '+str(j2)+' '+str(k2)+' xyz '+str(bounds[1][i1])+' '+str(bounds[1][j1])+' '+str(bounds[1][k1])+'\n')
                k1 = 5
                k2 = 5
            j1 = 4
            j2 = 5
        i1 = 3
        i2 = 4
    '''
    tg_lines.append('\nc Make block boundary constraints\n')
    blk_cnt = int((c[3]+1)/2)
    bbid = blk_cnt+(blk_cnt-1)*18
    tg_lines.append('bb 1 4 3 3 4 4 '+str(bbid)+';\n')
    tg_lines.append('bb 4 4 3 6 4 4 '+str(bbid+1)+';\n')
    tg_lines.append('bb 4 3 3 6 3 4 '+str(bbid+2)+';\n')
    tg_lines.append('bb 1 3 3 3 3 4 '+str(bbid+3)+';\n')
    tg_lines.append('bb 3 4 3 3 6 4 '+str(bbid)+';\n')
    tg_lines.append('bb 4 4 3 4 6 4 '+str(bbid+1)+';\n')
    tg_lines.append('bb 4 1 3 4 3 4 '+str(bbid+2)+';\n')
    tg_lines.append('bb 3 1 3 3 3 4 '+str(bbid+3)+';\n')
    tg_lines.append('bb 3 4 1 4 4 3 '+str(bbid+4)+';\n')
    tg_lines.append('bb 3 4 4 4 4 6 '+str(bbid+5)+';\n')
    tg_lines.append('bb 3 3 4 4 3 6 '+str(bbid+6)+';\n')
    tg_lines.append('bb 3 3 1 4 3 3 '+str(bbid+7)+';\n')
    tg_lines.append('bb 3 4 3 4 6 3 '+str(bbid+4)+';\n')
    tg_lines.append('bb 3 4 4 4 6 4 '+str(bbid+5)+';\n')
    tg_lines.append('bb 3 1 4 4 3 4 '+str(bbid+6)+';\n')
    tg_lines.append('bb 3 1 3 4 3 3 '+str(bbid+7)+';\n')
    tg_lines.append('bb 1 3 4 3 4 4 '+str(bbid+8)+';\n')
    tg_lines.append('bb 4 3 4 6 4 4 '+str(bbid+9)+';\n')
    tg_lines.append('bb 4 3 3 6 4 3 '+str(bbid+10)+';\n')
    tg_lines.append('bb 1 3 3 3 4 3 '+str(bbid+11)+';\n')
    tg_lines.append('bb 3 3 4 3 4 6 '+str(bbid+8)+';\n')
    tg_lines.append('bb 4 3 4 4 4 6 '+str(bbid+9)+';\n')
    tg_lines.append('bb 4 3 1 4 4 3 '+str(bbid+10)+';\n')
    tg_lines.append('bb 3 3 1 3 4 3 '+str(bbid+11)+';\n')
    tg_lines.append('bb 1 3 3 1 4 4 '+str(bbid+12)+';\n')
    tg_lines.append('bb 6 3 3 6 4 4 '+str(bbid+13)+';\n')
    tg_lines.append('bb 3 1 3 4 1 4 '+str(bbid+14)+';\n')
    tg_lines.append('bb 3 6 3 4 6 4 '+str(bbid+15)+';\n')
    tg_lines.append('bb 3 3 1 4 4 1 '+str(bbid+16)+';\n')
    tg_lines.append('bb 3 3 6 4 4 6 '+str(bbid+17)+';\n')
    
    tg_lines.append('\nc Surface Projections\n')
    tg_lines.append('sfi -1 0 -6;;;sd '+str(c[3]+1)+'\n')
    tg_lines.append('sfi ;; -1 0 -6;sd '+str(c[3]+1)+'\n')
    tg_lines.append('sfi ; -1 0 -6;;sd '+str(c[3]+1)+'\n')
    tg_lines.append('sfi -2 0 -5;;;sd '+str(c[3])+'\n')
    tg_lines.append('sfi ; -2 0 -5;;sd '+str(c[3])+'\n')
    tg_lines.append('sfi ;; -2 0 -5;sd '+str(c[3])+'\n')
    
    tg_lines.append('\nc Uniform smoothing\n')
    tg_lines.append('unifm 1 3 3 1 4 4 & 6 3 3 6 4 4 & 3 1 3 4 1 4 & 3 6 3 4 6 4 & 3 3 1 4 4 1 & 3 3 6 4 4 6 20 0 1;\n')
    tg_lines.append('unifm 1 3 3 6 4 4 & 3 4 3 4 6 4 & 3 1 3 4 3 4 & 3 3 1 4 4 3 & 3 3 4 4 4 6 20 0 1;\n')
            
    return tg_lines
    
     
    
        
    
if __name__ == '__main__':
    main(sys.argv[-1])