#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Created on Sat Jan  9 20:21:38 2016

@author: Edgar Haener
edgar.haner@gmail.com

Script to write poly file for Half-Cylinder device to a file that can be read by triangle.

"""
from __future__ import absolute_import, division, print_function
import matplotlib.pyplot as plt
import numpy as np


def writeFile():
    path='/home/magda/oomph-lib/trunk/user_drivers/edgar_haener/traingleMeshGen/'
    path='traingle/'
    
    baseUnit = 1
    w_channel= 4*baseUnit #Width main channel 
    r_obstacle = w_channel/4.0 #Radius of obstacle
    
    l_channel = 3.0 * w_channel + r_obstacle #length main channel
    w_diffuser = 6.75*baseUnit # width of the part of the diffuser angled at 45 degrees
    l_diffuser = 2 * baseUnit # final length
    nr_pts_obst = 100 # place a hundred points on the curved surface of the obstacle


    #Define the points on the boundary based on these measurments   
    A = np.array([0.0, -w_channel/2.0])
    B = A + np.array([0.0, w_channel])
    C = B + np.array([l_channel, 0.0])
    D = C + np.array([w_diffuser, w_diffuser])
    E = D + np.array([l_diffuser, 0.0 ])
    F = E +np.array([0.0, -2*w_diffuser-w_channel])
    G = F +np.array([-l_diffuser, 0.0])
    H = G + np.array([-w_diffuser, w_diffuser])
    
    #Obstacle
    I = H + np.array([0.0, w_channel/2.0 - r_obstacle])
    J = I + np.array([0.0, 2*r_obstacle])
    
    listP = [A, B, C, D, E, F, G, H, I, J]
    
    #now need to approximate the curved surface of the cylinder by a series of points
    angel_increment = np.pi / nr_pts_obst + 0.0
       
    theta=0.0
    for jj in range(nr_pts_obst-1):
        theta += angel_increment
        listP.append(np.array([l_channel - r_obstacle*np.sin(theta), r_obstacle*np.cos(theta)]))
    
    
    fig1 = plt.figure()
    for ii in range(8):
        if ii == 7:
            a=listP[ii]
            b=listP[0]
        else:
            a=listP[ii]
            b=listP[ii+1]
        plt.plot([a[0], b[0] ] , [a[1], b[1]], 'ko-')

    for ii in range(8, len(listP)):
        if ii == len(listP)-1:
            a=listP[ii]
            b=listP[8]
        else:
            a=listP[ii]
            b=listP[ii+1]
        plt.plot([a[0], b[0] ] , [a[1], b[1]], 'ko-')

    plt.show()    
    


    f = open(path + 'halfcylinder.poly', 'w')
    f.write('%d 2 0 0 # of pts, 2D, no attributes or boundary markers for points \n' %(len(listP)))
    
    #Writting points to file
    for ii in range(len(listP)):
        f.write('%d  %.4f %.4f \n' %(ii+1, listP[ii][0], listP[ii][1]))
    
    f.write('# END_OF_NODE_BLOCK \n\n')
    
    #writting connections
    
    f.write('%d 1  #   [number of segments (i.e. boundary edges), flag for using boundary markers here] \n' %(len(listP)))
    f.write('1 1 2  1 \n')
    f.write('2 2 3  2 \n')
    f.write('3 3 4  2 \n')
    f.write('4 4 5  2 \n')
    f.write('5 5 6  3 \n')
    f.write('6 6 7  4 \n')
    f.write('7 7 8  4 \n')
    f.write('8 8 1  4 \n')
    #Now obstacle
    f.write('9 9 10  5  \n')
    for kk in range(10, len(listP)):
        f.write('%d %d %d  5  \n' %(kk, kk, kk+1))
    f.write('%d %d %d  5  \n' %(len(listP), len(listP), 9))
    
    f.write('# END_OF_SEGMENT_BLOCK \n\n')
    
    #Holes
    f.write('1 #  [number of holes] \n')
    ypos = 0.0
    xpos= l_channel - r_obstacle/2.0
    f.write('1 %.4f %.4f \n' %(xpos, ypos))
    f.write(' # END_OF_HOLE_BLOCK\n')
    sV = '#baseUnit = %.2f , w_channel = %.2f ,  l_channel  = %.2f r_obstacle = %.2f , w_diffuser= %.2f, l_diffuser= %.2f \n' %(baseUnit, w_channel, l_channel, r_obstacle, w_diffuser, l_diffuser) 
    f.write('\n# Writen with Python Script: \n%s' %sV)
    
if __name__ =="__main__":
    writeFile()
#    plt.close('all')