Index: /anuga_core/source/anuga/pmesh/mesh_interface.py
===================================================================
--- /anuga_core/source/anuga/pmesh/mesh_interface.py	(revision 7698)
+++ /anuga_core/source/anuga/pmesh/mesh_interface.py	(revision 7699)
@@ -5,4 +5,5 @@
 import numpy as num
 from anuga.utilities.polygon import inside_polygon
+from anuga.utilities.polygon import polylist2points_verts
 import anuga.utilities.log as log
 
@@ -29,4 +30,5 @@
                              minimum_triangle_angle=28.0,
                              fail_if_polygons_outside=True,
+                             breaklines=None,
                              use_cache=False,
                              verbose=True):
@@ -65,4 +67,8 @@
     the lower left hand corner of  bounding_polygon (absolute)
     as the x and y values for the geo_ref.
+
+    breaklines is a list of polygons. These lines will be preserved by the
+               triangulation algorithm - useful for coastlines, walls, etc.
+               The polygons are not closed.			   
     
     Returns the mesh instance if no filename is given
@@ -90,4 +96,5 @@
               'minimum_triangle_angle': minimum_triangle_angle,
               'fail_if_polygons_outside': fail_if_polygons_outside,
+              'breaklines': breaklines,
               'verbose': verbose}   # FIXME (Ole): Should be bypassed one day. See ticket:14
 
@@ -124,4 +131,5 @@
                               minimum_triangle_angle=28.0,
                               fail_if_polygons_outside=True,
+                              breaklines=None,
                               verbose=True):
     """_create_mesh_from_regions - internal function.
@@ -216,5 +224,4 @@
 #    print 'check %s in production directory' %figname
 #    import sys; sys.exit()
-    
 
     if interior_holes is not None:        
@@ -246,4 +253,9 @@
 
     m = Mesh(geo_reference=mesh_geo_reference)
+
+    # build a list of discrete segments from the breakline polygons
+    if breaklines is not None:
+        points, verts = polylist2points_verts(breaklines)
+        m.add_points_and_segments(points, verts)
 
     # Do bounding polygon
@@ -303,5 +315,4 @@
 
 
-
     # NOTE (Ole): This was moved here as it is annoying if mesh is always
     # stored irrespective of whether the computation
Index: /anuga_core/source/anuga/pmesh/test_mesh_interface.py
===================================================================
--- /anuga_core/source/anuga/pmesh/test_mesh_interface.py	(revision 7698)
+++ /anuga_core/source/anuga/pmesh/test_mesh_interface.py	(revision 7699)
@@ -249,6 +249,6 @@
         self.failUnless(m.geo_reference.get_zone()==zone, 'FAILED!')
         self.failUnless(m.geo_reference.get_xllcorner()==min_x, 'FAILED!')
-        self.failUnless(m.geo_reference.get_yllcorner()==min_y, 'FAILED!')
-
+        self.failUnless(m.geo_reference.get_yllcorner()==min_y, 'FAILED!')		
+		
     def test_create_mesh_from_regions3(self):
         # These are the absolute values
@@ -636,7 +636,20 @@
             raise Exception, msg
             
-                    
-
-
+    def test_create_mesh_with_breaklines(self):
+         # These are the absolute values
+        polygon = [[100,100], [1000,100], [1000,1000], [100,1000]]
+
+        boundary_tags = {'walls': [0,1], 'bom': [2,3]}
+
+        m = create_mesh_from_regions(polygon,
+                                     boundary_tags,
+                                     10000000,
+                                     breaklines=[[[50,50],[2000,2000]]])
+									 
+        self.failUnless(len(m.regions) == 1, 'FAILED!')
+        segs = m.getUserSegments()
+        self.failUnless(len(segs) == 5, 'FAILED!')
+        self.failUnless(len(m.userVertices) == 6, 'FAILED!')
+		
     def test_create_mesh_from_regions_with_duplicate_verts(self):
         # These are the absolute values
Index: /anuga_core/source/anuga/utilities/polygon.py
===================================================================
--- /anuga_core/source/anuga/utilities/polygon.py	(revision 7698)
+++ /anuga_core/source/anuga/utilities/polygon.py	(revision 7699)
@@ -1355,6 +1355,21 @@
                               atol)
 
+
     return interpolated_values
 
+    
+def polylist2points_verts(polylist):
+    """ Convert a list of polygons to discrete points and vertices.
+    """
+    
+    offset = 0
+    points = []
+    vertices = []
+    for poly in polylist:
+        points.extend(poly)
+        vertices.extend([[i, i+1] for i in range(offset, offset+len(poly)-1)])
+        offset += len(poly)
+		
+    return points, vertices
 ##
 # @brief 
