Changeset 7699

Timestamp:

Apr 28, 2010, 5:32:21 PM (15 years ago)

Author:

James Hudson

Message:

Added breaklines to high level interface for task 305.

Location:

anuga_core/source/anuga

Files:

• pmesh/mesh_interface.py (modified) (8 diffs)
• pmesh/test_mesh_interface.py (modified) (2 diffs)
• utilities/polygon.py (modified) (1 diff)

anuga_core/source/anuga/pmesh/mesh_interface.py

@@ -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 +30 @@
                              minimum_triangle_angle=28.0,
                              fail_if_polygons_outside=True,
+                             breaklines=None,
                              use_cache=False,
                              verbose=True):
@@ -65 +67 @@
     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 +96 @@
               '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 +131 @@
                               minimum_triangle_angle=28.0,
                               fail_if_polygons_outside=True,
+                              breaklines=None,
                               verbose=True):
     """_create_mesh_from_regions - internal function.
@@ -216 +224 @@
 #    print 'check %s in production directory' %figname
 #    import sys; sys.exit()
-    
 
     if interior_holes is not None:        
@@ -246 +253 @@
 
     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 +315 @@
 
 
-
     # NOTE (Ole): This was moved here as it is annoying if mesh is always
     # stored irrespective of whether the computation

anuga_core/source/anuga/pmesh/test_mesh_interface.py

@@ -249 +249 @@
         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 +636 @@
             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

anuga_core/source/anuga/utilities/polygon.py

@@ -1355 +1355 @@
                               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