"""Class Domain - 1D domains for finite-volume computations of
   the shallow water wave equation


   Copyright 2004
   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
   Geoscience Australia
"""

class Domain:

    def __init__(self, coordinates, conserved_quantities = None,
                 other_quantities = None):
        """
        Build 1D elements from x coordinates
        """

        from Numeric import array, zeros, Float, Int

        self.beta = 1.0
        #Store Points
        self.coordinates = array(coordinates)

        #Register number of Elements
        self.number_of_elements = N = len(self.coordinates)-1
        

        #Allocate space for geometric quantities
        self.vertices  = zeros((N, 2), Float)
        self.centroids = zeros(N, Float)
        self.areas     = zeros(N, Float)
        for i in range(N):
            xl = self.coordinates[i]
            xr = self.coordinates[i+1]
            self.vertices[i,0] = xl
            self.vertices[i,1] = xr
            
            centroid = (xl+xr)/2
            self.centroids[i] = centroid

            msg = 'Coordinates should be ordered, smallest to largest'
            assert xr>xl, msg

            self.areas[i] = (xr-xl)

##         print 'N', N
##         print 'Centroid', self.centroids
##         print 'Areas', self.areas
##         print 'Vertex_Coordinates', self.vertices
        
    def get_centroids(self):
        """Return all coordinates of centroids
        Return x coordinate of centroid for each element as a N array
        """
        
        return self.centroids

    def get_vertices(self):
        """Return all coordinates of centroids
        Return x coordinate of centroid for each element as a N array
        """

        return self.vertices 

    def get_coordinate(self, elem_id, vertex=None):
        """Return coordinate of centroid,
        or left or right vertex. 
        Left vertex (vertex=0). Right vertex (vertex=1)
        """

        if vertex is None:
            return self.centroids[elem_id]
        else:
            return self.vertices[elem_id,vertex]

    def get_area(self, elem_id):
        """Return area of element id
        """

        return self.areas[elem_id]
    

if __name__ == "__main__":

    points1 = [0.0, 1.0, 2.0, 3.0]
    D1 = Domain(points1)

    print D1.get_coordinate(0)
    print D1.get_coordinate(0,1)

    try:
        print D1.get_coordinate(3)
    except:
        pass
    else:
        msg =  'Should have raised an out of bounds exception'
    raise msg
    
    #points2 = [0.0, 1.0, 2.0, 3.0, 2.5]
    #D2 = Domain(points2)