source: inundation/ga/storm_surge/pyvolution/generic_boundary_conditions.py @ 198

"""boundary.py - Classes for implementing boundary conditions
"""

from util import rotate

class Boundary:
    """Base class for boundary conditions.
       Specific boundary conditions must provide values for
       the conserved_quantities

       A boundary object has one neighbour; the one it
       serves.
    """   

    def __init__(self):
        pass
        
    def evaluate(self, vol_id=None, edge_id=None):
        msg = 'Generic class Boundary must be subclassed'
        raise msg


class Transmissive_boundary(Boundary):
    """Transmissive boundary returns same conserved quantities as
    those present in its neighbour volume.

    Underlying domain must be specified when boundary is instantiated
    """
    
    def __init__(self, domain = None):
        Boundary.__init__(self)

        if domain is None:
            msg = 'Domain must be specified for transmissive boundary'
            raise msg
        
        self.domain = domain
        
    def __repr__(self):
        return 'Transmissive_boundary(%s)' %self.domain

    def evaluate(self, vol_id, edge_id):
        """Transmissive boundaries return the edge values
        of the volume they serve.
        """
        
        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
        return q        


class Dirichlet_boundary(Boundary):
    """Dirichlet boundary returns constant values for the
    conserved quantities
    """


    def __init__(self, conserved_quantities=None):
        Boundary.__init__(self)

        if conserved_quantities is None:
            msg = 'Must specify one value for each conserved quantity'
            raise msg

        from Numeric import array, Float
        self.conserved_quantities=array(conserved_quantities).astype(Float)

    def __repr__(self):
        return 'Dirichlet boundary (%s)' %self.conserved_quantities
    
    def evaluate(self, vol_id=None, edge_id=None):
        return self.conserved_quantities



class Time_boundary(Boundary):
    """Time dependent boundary returns values for the
    conserved quantities as a function of time.
    Must specify domain to get access to model time and a function f(t)
    which must return conserved quantities as a function time
    """


    def __init__(self, domain = None, f=None):
        Boundary.__init__(self)

        try:
            q = f(0.0)
        except:
            raise 'Function for time boundary could not be executed'

        from Numeric import array, Float        
        try:
            q = array(q).astype(Float)
        except:
            msg = 'Return value from time boundary function could '
            msg += 'not be converted into a Numeric array of floats.\n'
            msg += 'Specified function should return either list or array.' 
            raise msg

        msg = 'ERROR: Time boundary function must return a 1d list or array '
        assert len(q.shape) == 1, msg
            
        d = len(domain.conserved_quantities)
        msg = 'Return value for f must be a list or an array of length %d' %d
        assert len(q) == d, msg

        self.f = f
        self.domain = domain

    def __repr__(self):
        return 'Time boundary'

    def evaluate(self, vol_id=None, edge_id=None):
        return self.f(self.domain.time)



#THIS FAR (10/8/4)
class Connective_boundary(Boundary):
    """Connective boundary returns values for the
    conserved quantities from a volume as defined by a connection table
    mapping between tuples of (volume id, face id) for volumes that
    have their boundaries connected. 

    FIXME: Perhaps include possibility of mapping between
    different domains as well

    FIXME: In case of shallow water we may want to have a
    special version that casts this in terms of height rather than stage
    """


    def __init__(self, table):
        from domain import Volume
        
        Boundary.__init__(self)

        self.connection_table = table
        self.Volume = Volume

    def __repr__(self):
        return 'Connective boundary'        

    #FIXME: IF we ever need to get field_values from connected volume,
    #that method could be overridden here (using same idea as in
    #get_conserved_quantities
    #def get_field_values()
    
    def get_conserved_quantities(self, volume, face=0):

        id = volume.id
        if self.connection_table.has_key((id, face)):
            other_id, other_face = self.connection_table[(id, face)] 
            
            other_volume = self.Volume.instances[other_id]
            cmd = 'q = other_volume.conserved_quantities_face%d' %face;
            exec(cmd)        
            return q
        else:
            msg = 'Volume, face tuple ($d, %d) has not been mapped'\
                  %(id, face)
            raise msg
            
        



#FIXME: Add a boundary with a general function of x,y, and t 

#FIXME: Add periodic boundaries e.g.:
# Attempt at periodic conditions from advection_spik. Remember this
#
#first = 2*(N-1)*N
#for i in range(1,2*N+1,2):
#    k = first + i-1#
#
#    print i,k
#    
#    domain[i].faces[2].neighbour = domain[k].faces[1]
#    domain[k].faces[1].neighbour = domain[i].faces[2]    
    


class General_boundary(Boundary):
    """General boundary which can compute conserved quantities based on
    their previous value, conserved quantities of its neighbour and model time.

    Must specify initial conserved quantities,
    neighbour,
    domain to get access to model time 
    a function f(q_old, neighbours_q, t) which must return
    new conserved quantities q as a function time

    FIXME: COMPLETE UNTESTED - JUST AN IDEA
    """

    def __init__(self, neighbour=None, conserved_quantities=None, domain=None, f=None):
        Boundary.__init__(self, neighbour=neighbour, conserved_quantities=conserved_quantities)

        self.f = f
        self.domain = domain


    def get_conserved_quantities(self, volume=None, face=0):
        return self.f(self.conserved_quantities,
                      neighbour.conserved_quantities,
                      self.domain.time)