Index: /trunk/anuga_core/source/anuga/operators/set_elevation.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/set_elevation.py	(revision 8931)
+++ /trunk/anuga_core/source/anuga/operators/set_elevation.py	(revision 8931)
@@ -0,0 +1,247 @@
+"""
+Set elevation operators
+
+
+"""
+
+__author__="steve"
+__date__ ="$09/03/2012 4:46:39 PM$"
+
+
+from anuga import Domain
+from anuga import Quantity
+import numpy as num
+import anuga.utilities.log as log
+
+from anuga.geometry.polygon import inside_polygon
+
+from anuga.operators.set_quantity import Set_quantity
+
+from anuga import indent
+
+class Set_elevation(Set_quantity):
+    """
+    Helper class to setup calculation of elevation
+    associated with a region (defined by indices, polygon or center/radius
+    """
+    
+    set_elevation = Set_quantity.set_value
+
+
+    def __init__(self,
+                 domain,
+                 elevation=None,
+                 indices=None,
+                 polygon=None,
+                 center=None,
+                 radius=None):
+
+
+        Set_quantity.__init__(self, domain, 'elevation',
+                              value = elevation,
+                              indices=indices,
+                              polygon=polygon,
+                              center=center,
+                              radius=radius,
+                              test_elevation=False)
+
+
+
+        #-----------------------------------------
+        # Extra structures to support maintaining
+        # continuity of elevation
+        #-----------------------------------------
+        self.setup_node_structures()
+
+
+        #------------------------------------------
+        # Extra aliases for changing elevation
+        # to ensure conservation of mass and
+        # continuity at vertices and edges
+        #------------------------------------------
+        self.elev_v  = self.domain.quantities['elevation'].vertex_values
+        self.elev_e = self.domain.quantities['elevation'].edge_values
+        self.stage_c = self.domain.quantities['stage'].centroid_values
+        self.elev_c = self.domain.quantities['elevation'].centroid_values
+
+        #------------------------------------------
+        # x,y coordinates of vertices of cells that are
+        # updated
+        #------------------------------------------
+        coord_v = self.domain.vertex_coordinates
+        ids = self.indices
+        if ids is None:
+            self.v_x = coord_v[:,0].reshape((-1,3))
+            self.v_y = coord_v[:,1].reshape((-1,3))
+        else:
+            self.v_x = coord_v[:,0].reshape((-1,3))[ids]
+            self.v_y = coord_v[:,1].reshape((-1,3))[ids]
+
+        #print self.v_x.shape
+        #print self.v_y.shape
+
+        #------------------------------------------
+        # Need to turn off this optimization as it
+        # doesn't fixup the relationship between
+        # bed and stage vertex values in dry region
+        #------------------------------------------
+        self.domain.optimise_dry_cells = 0
+
+
+
+
+    def __call__(self):
+        """
+        Apply rate to those triangles defined in indices
+
+        indices == [], then don't apply anywhere
+        indices == None, then apply everywhere
+        otherwise apply for the specific indices
+        """
+
+        if self.value is None:
+            return
+
+        if self.indices is []:
+            return
+
+        #------------------------------------------
+        # Apply changes to elevation vertex values
+        # via the update_quantites routine
+        #------------------------------------------
+        if not self.update_quantities():
+            return
+
+
+        #------------------------------------------
+        # Cleanup elevation and stage quantity values
+        #-----------------------------------------
+        if self.indices is None:
+
+            #--------------------------------------
+            # Make elevation continuous and clean up
+            # stage values to ensure conservation
+            #--------------------------------------
+            height_c = self.stage_c - self.elev_c
+            self.domain.quantities['elevation'].smooth_vertex_values()
+            self.domain.quantities['elevation'].interpolate()
+            self.stage_c[:] = self.elev_c +  height_c
+
+
+        else:
+
+            #--------------------------------------
+            # Make elevation continuous and clean up
+            # stage values to ensure conservation
+            #--------------------------------------
+            height_c = self.stage_c[self.vols] - self.elev_c[self.vols]
+            for nid in self.node_ids:
+                non = self.domain.number_of_triangles_per_node[nid]
+
+                vid = num.arange(self.node_index[nid], self.node_index[nid+1],dtype=num.int)
+                vidd = self.domain.vertex_value_indices[vid]
+
+                self.elev_v[vidd/3,vidd%3] = num.sum(self.elev_v[vidd/3,vidd%3])/non
+
+
+            #--------------------------------------
+            # clean up the centroid values and edge values
+            #--------------------------------------
+            self.elev_c[self.vols] = num.mean(self.elev_v[self.vols],axis=1)
+
+            self.elev_e[self.vols,0] = 0.5*(self.elev_v[self.vols,1]+ self.elev_v[self.vols,2])
+            self.elev_e[self.vols,1] = 0.5*(self.elev_v[self.vols,2]+ self.elev_v[self.vols,0])
+            self.elev_e[self.vols,2] = 0.5*(self.elev_v[self.vols,0]+ self.elev_v[self.vols,1])
+
+            self.stage_c[self.vols] = self.elev_c[self.vols] +  height_c
+
+
+
+    def update_quantities(self):
+        """Update the vertex values of the quantities to model erosion
+        """
+
+        if self.value is None:
+            return False
+
+
+        updated = True
+
+        if self.indices is None:
+
+            #--------------------------------------
+            # Update all three vertices for each cell
+            #--------------------------------------
+            self.elev_v[:] = self.get_value(self.v_x, self.v_y)
+
+        else:
+
+            #--------------------------------------
+            # Update all three vertices for each cell
+            # cell associated with ids
+            #--------------------------------------
+            ids = self.indices
+            self.elev_v[ids] = self.get_value(self.v_x, self.v_y)
+
+
+        return updated
+
+
+
+
+    def setup_node_structures(self):
+        """ For setting elevation we need to
+        ensure that the elevation quantity remains
+        continuous (at least for version 1.3 of anuga)
+
+        So we need to find all the vertices that need to
+        update within each timestep.
+        """
+
+
+
+        if self.indices is None or self.indices is []:
+            self.vol_ids  = None
+            self.vols = None
+            self.vert_ids = None
+            return
+
+
+        
+        node_ids = set()
+
+        for ind in self.indices:
+            for k in [0,1,2]:
+                node_ids.add(self.domain.triangles[ind,k])
+
+        self.node_ids = [ id for id in node_ids ]
+
+
+
+        node_index = num.zeros((self.domain.number_of_nodes)+1, dtype = num.int)
+
+        # FIXME: SR Don't we calculate this for the domain already!
+        k = 0
+        node_index[0] = 0
+        for i in range(self.domain.number_of_nodes):
+            node_index[i+1] = node_index[i] + self.domain.number_of_triangles_per_node[i]
+
+        self.node_index = node_index
+
+        vertex_ids =[]
+        for nid in self.node_ids:
+            #print nid,self.domain.number_of_triangles_per_node[nid]
+            for vid in range(node_index[nid], node_index[nid+1]):
+                vidd = self.domain.vertex_value_indices[vid]
+                vertex_ids.append(vidd)
+                #print '   ',nid, vid, vidd, vidd/3, vidd%3
+
+        self.vol_ids  = num.array(vertex_ids,dtype=num.int)/3
+        self.vols = num.array(list(set(self.vol_ids)), dtype=num.int)
+        self.vert_ids = num.array(vertex_ids,dtype=num.int)%3
+
+
+
+
+
+
Index: /trunk/anuga_core/source/anuga/operators/set_quantity.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/set_quantity.py	(revision 8931)
+++ /trunk/anuga_core/source/anuga/operators/set_quantity.py	(revision 8931)
@@ -0,0 +1,148 @@
+"""
+Set elevation operators
+
+
+"""
+
+__author__="steve"
+__date__ ="$09/03/2012 4:46:39 PM$"
+
+import numpy as num
+
+from anuga import Domain
+from anuga import Quantity
+
+import anuga.utilities.log as log
+from anuga.geometry.polygon import inside_polygon
+from anuga.utilities.function_utils import determine_function_type
+from anuga.operators.region import Region
+from anuga import indent
+
+class Set_quantity(Region):
+    """
+    Helper class to setup calculation of quantity
+    associated with a region (defined by indices, polygon or center/radius
+    """
+
+    def __init__(self,
+                 domain,
+                 quantity,
+                 value=None,
+                 indices=None,
+                 polygon=None,
+                 center=None,
+                 radius=None,
+                 verbose = False,
+                 test_elevation=True):
+
+
+        Region.__init__(self, domain,
+                        indices=indices,
+                        polygon=polygon,
+                        center=center,
+                        radius=radius,
+                        verbose=verbose)
+        
+
+        self.set_value(value)
+
+        #-------------------------------------------
+        # Test quantity
+        #-------------------------------------------
+        self.quantity = quantity
+        msg = 'quantity not found in domain'
+        assert quantity in domain.quantities, msg
+
+        if test_elevation:
+            msg ='Use Set_elevation to maintain continuity'
+            assert quantity is not 'elevation'
+        
+        #-------------------------------------------
+        # Useful quantity alias
+        #------------------------------------------
+        self.quantity_c = self.domain.quantities[quantity].centroid_values
+
+        self.coord_c = self.domain.centroid_coordinates
+        self.areas = self.domain.areas
+
+
+    def __call__(self):
+        """
+        Apply value to those triangles defined by indices
+
+        indices == [], don't apply anywhere
+        indices == None, apply everywhere
+        otherwise apply for the specific indices
+        """
+
+        if self.indices is []:
+            return
+
+        value = self.get_value()
+
+        if value is None:
+            return
+
+
+        if self.indices is None:
+
+            #--------------------------------------
+            # Update all three vertices for each cell
+            #--------------------------------------
+            self.quantity_c[:] = self.get_value(self.coord_c[:,0], self.coord_c[:,1])
+
+        else:
+
+            #--------------------------------------
+            # Update all three vertices for each cell
+            #--------------------------------------
+            ids = self.indices
+            x = self.coord_c[ids,0]
+            y = self.coord_c[ids,1]
+            self.quantity_c[ids] = self.get_value(x,y)
+
+
+
+    def set_value(self, value = None):
+
+        self.value = value
+        self.value_type = determine_function_type(value)
+
+
+        
+    def get_value(self, x = None, y = None, t = None):
+        """Get value of quantity at time t.
+        If t not specified, return value at current domain time
+        """
+
+        #from anuga.fit_interpolate.interpolate import Modeltime_too_early, \
+        #                                              Modeltime_too_late
+
+
+        if t is None:
+            t = self.domain.get_time()
+
+        #try:
+        if self.value_type == 't':
+            value = self.value(t)
+        elif self.value_type == 'x,y':
+            value = self.value(x,y)
+        elif self.value_type == 'x,y,t':
+            value = self.value(x,y,t)
+        else:
+            value = float(self.value)
+
+        #except Modeltime_too_early, e:
+        #    raise Modeltime_too_early(e)
+        #except Modeltime_too_late, e:
+        #    raise Modeltime_too_late(e)
+
+
+        return value
+
+
+
+
+
+
+
Index: /trunk/anuga_core/source/anuga/operators/set_quantity_operator.py
===================================================================
--- /trunk/anuga_core/source/anuga/operators/set_quantity_operator.py	(revision 8931)
+++ /trunk/anuga_core/source/anuga/operators/set_quantity_operator.py	(revision 8931)
@@ -0,0 +1,151 @@
+"""
+Set quantity operators
+
+
+"""
+
+__author__="steve"
+__date__ ="$09/03/2012 4:46:39 PM$"
+
+
+from anuga import Domain
+from anuga import Quantity
+import numpy as num
+import anuga.utilities.log as log
+
+from anuga.geometry.polygon import inside_polygon
+
+from anuga.operators.base_operator import Operator
+from anuga.operators.set_quantity import Set_quantity
+
+from anuga import indent
+
+
+#===============================================================================
+# General Set Quantity  Operator
+#===============================================================================
+
+class Set_quantity_operator(Operator, Set_quantity):
+    """
+    Set the elevation in a region (careful to maintain continuitiy of elevation)
+
+    indices: None == all triangles, Empty list [] no triangles
+
+    elevation can be a function of time.
+
+    """
+
+
+    __call__ = Set_quantity.__call__
+
+    
+    def __init__(self,
+                 domain,
+                 quantity,
+                 value=None,
+                 indices=None,
+                 polygon=None,
+                 center=None,
+                 radius=None,
+                 description = None,
+                 label = None,
+                 logging = False,
+                 verbose = False):
+
+
+ 
+
+        Set_quantity.__init__(self, domain, quantity, value, indices, polygon, center, radius)
+
+        Operator.__init__(self, domain, description, label, logging, verbose)
+
+
+
+
+    def parallel_safe(self):
+        """Operator is applied independently on each cell and
+        so is parallel safe.
+        """
+        return True
+
+    def statistics(self):
+
+        message = self.label + ': Set_elevation_operator'
+        message = message + ' on triangles '+ str(self.indices)
+        return message
+
+
+    def timestepping_statistics(self):
+
+        #message  = indent + self.label + ': Set_elevation = ' + str('')
+        #message  += ' at center '+str(self.center)
+
+        return 'Set_quantity_operator active at time '+str(self.domain.get_time())
+
+
+
+#===============================================================================
+# Specific Bed Operators for circular region.
+#===============================================================================
+class Circular_set_quantity_operator(Set_quantity_operator):
+    """
+    Set elevation over a circular region
+
+    """
+
+    def __init__(self, domain, quantity,
+                 value=0.0,
+                 center=None,
+                 radius=None,
+                 verbose=False):
+
+
+        assert center is not None
+        assert radius is not None
+
+
+        Set_quantity_operator.__init__(self,
+                                       domain,
+                                       quantity,
+                                       value=value,
+                                       center=center,
+                                       radius=radius,
+                                       verbose=verbose)
+
+
+
+
+
+#===============================================================================
+# Specific Bed Operators for polygonal region.
+#===============================================================================
+class Polygonal_set_quantity_operator(Set_quantity_operator):
+    """
+    Add water at certain rate (ms^{-1} = vol/Area/sec) over a
+    polygonal region
+
+    rate can be a function of time.
+
+    """
+
+    def __init__(self, domain, quantity,
+                 value=0.0,
+                 polygon=None,
+                 verbose=False):
+
+
+
+        Set_elevation_operator.__init__(self,
+                               domain,
+                               quantity,
+                               value=value,
+                               polygon=polygon,
+                               verbose=verbose)
+
+
+
+        
+
+
+
+
