source: branches/numpy/anuga/abstract_2d_finite_volumes/mesh_factory.py @ 6883

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

import numpy as num


def rectangular_old(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 anuga.config import epsilon

    deltax = float(len1)/m
    deltay = 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 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 anuga.config import epsilon

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

    #Calculate number of points
    Np = (m+1)*(n+1)

    class Index:

        def __init__(self, n,m):
            self.n = n
            self.m = m

        def __call__(self, i,j):
            return j+i*(self.n+1)


    index = Index(n,m)

    points = num.zeros((Np, 2), num.float)

    for i in range(m+1):
        for j in range(n+1):

            points[index(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]

    #Construct 2 triangles per rectangular element and assign tags to boundary
    #Calculate number of triangles
    Nt = 2*m*n


    elements = num.zeros((Nt, 3), num.int)
    boundary = {}
    nt = -1
    for i in range(m):
        for j in range(n):
            nt = nt + 1
            i1 = index(i,j+1)
            i2 = index(i,j)
            i3 = index(i+1,j+1)
            i4 = index(i+1,j)


            #Update boundary dictionary and create elements
            if i == m-1:
                boundary[nt, 2] = 'right'
            if j == 0:
                boundary[nt, 1] = 'bottom'
            elements[nt,:] = [i4,i3,i2] #Lower element
            nt = nt + 1

            if i == 0:
                boundary[nt, 2] = 'left'
            if j == n-1:
                boundary[nt, 1] = 'top'
            elements[nt,:] = [i1,i2,i3] #Upper element

    return points, elements, boundary


def rectangular_cross(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 anuga.config import epsilon

    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([delta1*i + origin[0], delta2*j  + origin[1]])

    # Construct 4 triangles per element
    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]
            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25

            # Create centre point
            v5 = len(points)
            points.append([x, y])

            #Create left triangle
            if i == 0:
                boundary[(len(elements), 1)] = 'left'
            elements.append([v2,v5,v1])

            #Create bottom triangle
            if j == 0:
                boundary[(len(elements), 1)] = 'bottom'
            elements.append([v4,v5,v2])

            #Create right triangle
            if i == m-1:
                boundary[(len(elements), 1)] = 'right'
            elements.append([v3,v5,v4])

            #Create top triangle
            if j == n-1:
                boundary[(len(elements), 1)] = 'top'
            elements.append([v1,v5,v3])


    return points, elements, boundary

def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
    """Setup a oblique grid of triangles
    with m segments in the x-direction and n segments in the y-direction

    """

    import math

    from anuga.config import epsilon


    deltax = lenx/float(m)
    deltay = leny/float(n)
    a = 0.75*lenx*math.tan(theta/180.*math.pi)
    x1 = lenx
    y1 = 0
    x2 = lenx
    y2 = leny
    x3 = 0.25*lenx
    y3 = leny
    x4 = x3
    y4 = 0
    a2 = a/(x1-x4)
    a1 = -a2*x4
    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3

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

    for i in range(m+1):
        x = deltax*i
        for j in range(n+1):
            y = deltay*j
            if x > 0.25*lenx:
                y = a1 + a2*x + a3*y + a4*x*y

            vertices[i,j] = len(points)
            points.append([x + origin[0], y + origin[1]])

    # Construct 2 triangles per element
    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

            if i == 0:
                boundary[(len(elements), 2)] = 'left'
            if j == n-1:
                boundary[(len(elements), 1)] = 'top'

            elements.append([v1,v2,v3]) #Upper

    return points, elements, boundary


def circular(m, n, radius=1.0, center = (0.0, 0.0)):
    """Setup a circular grid of triangles with m concentric circles and
    with n radial segments. If radius is are omitted the mesh defaults to
    the unit circle radius.

    radius: radius of circle

    #FIXME: The triangles become degenerate for large values of m or n.
    """



    from math import pi, cos, sin

    radius = float(radius) #Ensure floating point format

    #Dictionary of vertex objects and list of points
    vertices = {}
    points = [[0.0, 0.0]] #Center point
    vertices[0, 0] = 0

    for i in range(n):
        theta = 2*i*pi/n
        x = cos(theta)
        y = sin(theta)
        for j in range(1,m+1):
            delta = j*radius/m
            vertices[i,j] = len(points)
            points.append([delta*x, delta*y])

    #Construct 2 triangles per element
    elements = []
    for i in range(n):
        for j in range(1,m):

            i1 = (i + 1) % n  #Wrap around

            v1 = vertices[i,j+1]
            v2 = vertices[i,j]
            v3 = vertices[i1,j+1]
            v4 = vertices[i1,j]

            elements.append([v4,v2,v3]) #Lower
            elements.append([v1,v3,v2]) #Upper


    #Do the center
    v1 = vertices[0,0]
    for i in range(n):
        i1 = (i + 1) % n  #Wrap around
        v2 = vertices[i,1]
        v3 = vertices[i1,1]

        elements.append([v1,v2,v3]) #center

    return points, elements


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

    from anuga.utilities.numerical_tools 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, values



def strang_mesh(filename):
    """Read Strang generated mesh.
    """

    from math import pi
    from anuga.utilities.numerical_tools import anglediff


    fid = open(filename)
    points = []    # List of x, y coordinates
    triangles = [] # List of vertex ids as listed in the file

    for line in fid.readlines():
        fields = line.split()
        if len(fields) == 2:
            # we are reading vertex coordinates
            points.append([float(fields[0]), float(fields[1])])
        elif len(fields) == 3:
            # we are reading triangle point id's (format ae+b)
            triangles.append([int(float(fields[0]))-1,
                              int(float(fields[1]))-1,
                              int(float(fields[2]))-1])
        else:
            raise 'wrong format in ' + filename

    elements = [] #Final list of elements

    for t in triangles:
        #Get vertex coordinates
        v0 = t[0]
        v1 = t[1]
        v2 = t[2]

        x0 = points[v0][0]
        y0 = points[v0][1]
        x1 = points[v1][0]
        y1 = points[v1][1]
        x2 = points[v2][0]
        y2 = points[v2][1]

        #Check that points are arranged in counter clock-wise order
        vec0 = [x1-x0, y1-y0]
        vec1 = [x2-x1, y2-y1]
        vec2 = [x0-x2, y0-y2]

        a0 = anglediff(vec1, vec0)
        a1 = anglediff(vec2, vec1)
        a2 = anglediff(vec0, vec2)

        if a0 < pi and a1 < pi and a2 < pi:
            elements.append([v0, v1, v2])
        else:
            elements.append([v0, v2, v1])

    return points, elements

# #Map from edge number to indices of associated vertices
# edge_map = ((1,2), (0,2), (0,1))

def contracting_channel(m, n, W_upstream = 1., W_downstream = 0.75,
                            L_1 = 5.0, L_2 = 2.0, L_3 = 10, origin = (0.0, 0.0)):
    """Setup a contracting channel grid of triangles
    with m segments in the x-direction and n segments in the y-direction

    """

    import math

    from anuga.config import epsilon


    lenx = L_1 + L_2 + L_3
    leny = W_upstream
    deltax = lenx/float(m)
    deltay = leny/float(n)

    x1 = 0
    y1 = 0
    x2 = L_1
    y2 = 0
    x3 = L_1 + L_2
    y3 = (W_upstream - W_downstream)/2
    x4 = L_1 + L_2 + L_3
    y4 = y3
    x5 = x4
    y5 = y4 + W_downstream
    x6 = L_1 + L_2
    y6 = y5
    x7 = L_1
    y7 = W_upstream
    x8 = 0
    y8 = W_upstream
    a1 = 0
    a2 = (W_upstream - W_downstream)/(2*L_2)
    a3 = 1
    a4 = (W_downstream - W_upstream)/(L_2*W_upstream)

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

    for i in range(m+1):
        x = deltax*i
        for j in range(n+1):
            y = deltay*j
            if x > L_1 and x <= (L_1 + L_2):
                y = a1 + a2*(x - L_1) + a3*y + a4*(x - L_1)*y
            elif x > L_1 + L_2:
                y = (W_upstream - W_downstream)/2 + deltay*j*W_downstream/W_upstream

            vertices[i,j] = len(points)
            points.append([x + origin[0], y + origin[1]])

    # Construct 2 triangles per element
    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

            if i == 0:
                boundary[(len(elements), 2)] = 'left'
            if j == n-1:
                boundary[(len(elements), 1)] = 'top'

            elements.append([v1,v2,v3]) #Upper

    return points, elements, boundary


def contracting_channel_cross(m, n, W_upstream = 1., W_downstream = 0.75,
                              L_1 = 5.0, L_2 = 2.0, L_3 = 10, origin = (0.0, 0.0)):
    """Setup a contracting channel grid of triangles
    with m segments in the x-direction and n segments in the y-direction

    """

    import math

    from anuga.config import epsilon


    lenx = L_1 + L_2 + L_3
    leny = W_upstream
    deltax = lenx/float(m)
    deltay = leny/float(n)

    x1 = 0
    y1 = 0
    x2 = L_1
    y2 = 0
    x3 = L_1 + L_2
    y3 = (W_upstream - W_downstream)/2
    x4 = L_1 + L_2 + L_3
    y4 = y3
    x5 = x4
    y5 = y4 + W_downstream
    x6 = L_1 + L_2
    y6 = y5
    x7 = L_1
    y7 = W_upstream
    x8 = 0
    y8 = W_upstream
    a1 = 0
    a2 = (W_upstream - W_downstream)/(2*L_2)
    a3 = 1
    a4 = (W_downstream - W_upstream)/(L_2*W_upstream)

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

    for i in range(m+1):
        x = deltax*i
        for j in range(n+1):
            y = deltay*j
            if x > L_1 and x <= (L_1 + L_2):
                y = a1 + a2*(x - L_1) + a3*y + a4*(x - L_1)*y
            elif x > L_1 + L_2:
                y = (W_upstream - W_downstream)/2 + deltay*j*W_downstream/W_upstream

            vertices[i,j] = len(points)
            points.append([x + origin[0], y + origin[1]])

    # Construct 4 triangles per element
    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]
            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
            v5 = len(points)
            points.append([x, y])

            #Create left triangle
            if i == 0:
                boundary[(len(elements), 1)] = 'left'
            elements.append([v2,v5,v1])

            #Create bottom triangle
            if j == 0:
                boundary[(len(elements), 1)] = 'bottom'
            elements.append([v4,v5,v2])

            #Create right triangle
            if i == m-1:
                boundary[(len(elements), 1)] = 'right'
            elements.append([v3,v5,v4])

            #Create top triangle
            if j == n-1:
                boundary[(len(elements), 1)] = 'top'
            elements.append([v1,v5,v3])


    return points, elements, boundary




def oblique_cross(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
    """Setup a oblique grid of triangles
    with m segments in the x-direction and n segments in the y-direction

    """

    import math

    from anuga.config import epsilon


    deltax = lenx/float(m)
    deltay = leny/float(n)
    a = 0.75*lenx*math.tan(theta/180.*math.pi)
    x1 = lenx
    y1 = 0
    x2 = lenx
    y2 = leny
    x3 = 0.25*lenx
    y3 = leny
    x4 = x3
    y4 = 0
    a2 = a/(x1-x4)
    a1 = -a2*x4
    a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3)
    a3 = 1. - (a1 + a2*x3)/y3 - a4*x3

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

    for i in range(m+1):
        x = deltax*i
        for j in range(n+1):
            y = deltay*j
            if x > 0.25*lenx:
                y = a1 + a2*x + a3*y + a4*x*y

            vertices[i,j] = len(points)
            points.append([x + origin[0], y + origin[1]])

    # Construct 4 triangles per element
    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]
            x = (points[v1][0]+points[v2][0]+points[v3][0]+points[v4][0])*0.25
            y = (points[v1][1]+points[v2][1]+points[v3][1]+points[v4][1])*0.25
            v5 = len(points)
            points.append([x, y])

            #Update boundary dictionary and create elements
                        #Create left triangle
            if i == 0:
                boundary[(len(elements), 1)] = 'left'
            elements.append([v2,v5,v1])

            #Create bottom triangle
            if j == 0:
                boundary[(len(elements), 1)] = 'bottom'
            elements.append([v4,v5,v2])

            #Create right triangle
            if i == m-1:
                boundary[(len(elements), 1)] = 'right'
            elements.append([v3,v5,v4])

            #Create top triangle
            if j == n-1:
                boundary[(len(elements), 1)] = 'top'
            elements.append([v1,v5,v3])


    return points, elements, boundary