"""Library of standard meshes and facilities for reading various
mesh file formats
"""


def rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):

    """Setup a rectangular grid of triangles
    with m+1 by n+1 grid points
    and side lengths len1, len2. If side lengths are omitted
    the mesh defaults to the unit square.
    
    len1: x direction (left to right)
    len2: y direction (bottom to top)

    Return to lists: points and elements suitable for creating a Mesh or
    FVMesh object, e.g. Mesh(points, elements) 
    """

    from config import epsilon
    
    #E = m*n*2        #Number of triangular elements
    #P = (m+1)*(n+1)  #Number of initial vertices

    delta1 = float(len1)/m
    delta2 = float(len2)/n    

    #Dictionary of vertex objects
    vertices = {}
    points = []

    for i in range(m+1):
        for j in range(n+1):
            vertices[i,j] = len(points)
            points.append([i*delta1 + origin[0], j*delta2 + origin[1]])



    #Construct 2 triangles per rectangular element and assign tags to boundary
    elements = []
    boundary = {}
    for i in range(m):
        for j in range(n):
            v1 = vertices[i,j+1]
            v2 = vertices[i,j]            
            v3 = vertices[i+1,j+1]            
            v4 = vertices[i+1,j]            

            #Update boundary dictionary and create elements
            if i == m-1:
                boundary[(len(elements), 2)] = 'right'
            if j == 0:
                boundary[(len(elements), 1)] = 'bottom'                
            elements.append([v4,v3,v2]) #Lower element

            if i == 0:
                boundary[(len(elements), 2)] = 'left'
            if j == n-1:
                boundary[(len(elements), 1)] = 'top' 
            elements.append([v1,v2,v3]) #Upper element   

    return points, elements, boundary