def quad4_subdiv(elem_coords, intervals):
    """Compute an equispaced subdivision of a quad4 element.

    Quadrature points are equispaced, rather than Gaussian; improved
    order of accuracy of Gaussian quadrature is not achieved for
    discontinuous data.

    Input:
    elem_coords:  coordinates of element's nodes, assuming exodus
                    node order convention - (counter clockwise around
                    the element)
    intervals:    The element will be subdivided into nx*ny equispaced
                    subelements, where nx = ny = intervals
    """

    subelem_coords = []

    for jj in range(intervals):
        j = (0.5 + jj)/intervals
        for ii in range(intervals):
            i = (0.5 + ii)/intervals
            X_pt = (   (1-j)*(1-i)*elem_coords[0][0]
                     + (1-j)*   i *elem_coords[1][0]
                     +    j *   i *elem_coords[2][0]
                     +    j *(1-i)*elem_coords[3][0] )
            Y_pt = (   (1-j)*(1-i)*elem_coords[0][1]
                     + (1-j)*   i *elem_coords[1][1]
                     +    j *   i *elem_coords[2][1]
                     +    j *(1-i)*elem_coords[3][1] )
            subelem_coords.append( [X_pt, Y_pt, 0.0] )

    return subelem_coords