def triangle_quad(triangle, n=6):
    """
    Given a triangle in the form [(r1,z1), (r2,z2), (r3,z3)]
    returns a list of evaluation points [(r,z,weight),...]

    n   number of quadrature points.
        currently: 1, 3 or 6

    Coefficients taken from http://www.cs.rpi.edu/~flaherje/pdf/fea6.pdf
    Joseph E. Flaherty course notes, Rensselaer Polytechnic Institute
    """
    assert len(triangle) == 3
    r1, z1 = triangle[0]
    r2, z2 = triangle[1]
    r3, z3 = triangle[2]
    if n == 1:
        return [((r1 + r2 + r3) / 3, (z1 + z2 + z3) / 3, 1.0)]
    elif n == 3:
        return [((4 * r1 + r2 + r3) / 6, (4 * z1 + z2 + z3) / 6, 1.0 / 3),
            ((r1 + 4 * r2 + r3) / 6, (z1 + 4 * z2 + z3) / 6, 1.0 / 3), ((r1 +
            r2 + 4 * r3) / 6, (z1 + z2 + 4 * z3) / 6, 1.0 / 3)]
    elif n == 6:
        a = 0.816847572980459
        b = 0.5 * (1.0 - a)
        c = 0.10810301816807
        d = 0.5 * (1.0 - c)
        return [(a * r1 + b * r2 + b * r3, a * z1 + b * z2 + b * z3, 
            0.109951743655322), (b * r1 + a * r2 + b * r3, b * z1 + a * z2 +
            b * z3, 0.109951743655322), (b * r1 + b * r2 + a * r3, b * z1 +
            b * z2 + a * z3, 0.109951743655322), (c * r1 + d * r2 + d * r3,
            c * z1 + d * z2 + d * z3, 0.223381589678011), (d * r1 + c * r2 +
            d * r3, d * z1 + c * z2 + d * z3, 0.223381589678011), (d * r1 +
            d * r2 + c * r3, d * z1 + d * z2 + c * z3, 0.223381589678011)]
    else:
        raise ValueError('Quadrature not available for n={}'.format(n))
