source: storm_surge/pyvolution/mesh_factory.py @ 907

"""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       



def from_polyfile(name):
    """Read mesh from .poly file, an obj like file format
    listing first vertex coordinates and then connectivity
    """

    from util import anglediff
    from math import pi
    import os.path
    root, ext = os.path.splitext(name)

    if ext == 'poly':
        filename = name
    else:
        filename = name + '.poly'

    
    fid = open(filename)

    points = []    #x, y
    values = []    #z
    vertex_values = []    #Repeated z    
    triangles = [] #v0, v1, v2
    
    lines = fid.readlines()

    keyword = lines[0].strip()
    msg = 'First line in .poly file must contain the keyword: POINTS'
    assert keyword == 'POINTS', msg

    offending = 0
    i = 1
    while keyword == 'POINTS':
        line = lines[i].strip()
        i += 1

        if line == 'POLYS':
            keyword = line
            break

        fields = line.split(':')
        assert int(fields[0]) == i-1, 'Point indices not consecutive'

        #Split the three floats
        xyz = fields[1].split()

        x = float(xyz[0])
        y = float(xyz[1])
        z = float(xyz[2])        

        points.append([x, y])
        values.append(z)
        

    k = i    
    while keyword == 'POLYS':
        line = lines[i].strip()
        i += 1

        if line == 'END':
            keyword = line
            break
        

        fields = line.split(':')
        assert int(fields[0]) == i-k, 'Poly indices not consecutive'

        #Split the three indices
        vvv = fields[1].split()

        i0 = int(vvv[0])-1
        i1 = int(vvv[1])-1
        i2 = int(vvv[2])-1

        #Check for and exclude degenerate areas
        x0 = points[i0][0]
        y0 = points[i0][1]
        x1 = points[i1][0]
        y1 = points[i1][1]
        x2 = points[i2][0]
        y2 = points[i2][1]                

        area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
        if area > 0:

            #Ensure that points are arranged in counter clock-wise order
            v0 = [x1-x0, y1-y0]
            v1 = [x2-x1, y2-y1]
            v2 = [x0-x2, y0-y2]

            a0 = anglediff(v1, v0)
            a1 = anglediff(v2, v1)
            a2 = anglediff(v0, v2)        

            
            if a0 < pi and a1 < pi and a2 < pi:
                #all is well
                j0 = i0
                j1 = i1
                j2 = i2
            else:
                #Swap two vertices
                j0 = i1
                j1 = i0
                j2 = i2
                
            triangles.append([j0, j1, j2])
            vertex_values.append([values[j0], values[j1], values[j2]])
        else:
            offending +=1 

    print 'Removed %d offending triangles out of %d' %(offending, len(lines))
    return points, triangles, vertex_values