Changeset 3719

Timestamp:

Oct 10, 2006, 10:48:31 AM (17 years ago)

Author:

steve

Message:

Adding holes to create_mesh_from_regions

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity_ext.c (modified) (1 diff)
• examples/island.py (modified) (4 diffs)
• mesh_engine/test_triangmodule.py (modified) (1 diff)
• pmesh/mesh_interface.py (modified) (6 diffs)
• pmesh/test_mesh_interface.py (modified) (2 diffs)
• shallow_water/shallow_water_ext.c (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -209 +209 @@
 
 
-        //MH080605 set semi_impliit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
+        //MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
         for (k=0;k<N;k++){
                 semi_implicit_update[k]=0.0;

anuga_core/source/anuga/examples/island.py

@@ -33 +33 @@
                                            'top': [2],
                                            'left': [3]},
-                          maximum_triangle_area = 5,
+                          maximum_triangle_area = 1.0,
                           filename = 'island.msh' ,
-                          interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 3)]
+                          interior_regions=[ ([[50,25], [70,25], [70,75], [50,75]], 1.0)]
+                          #interior_holes=[[50,25], [70,25], [70,75], [50,75]],
                           )
 
@@ -54 +55 @@
 # beta_h == 0.0 means that constant (1st order) gradients are introduced
 # on h. This is equivalent to the constant depth used previously.
-domain.beta_h = 0.9
+domain.beta_h     = 0.2
+domain.beta_w_dry = 0.2
 
 
@@ -78 +80 @@
 
 #domain.set_quantity('friction', 0.1)  #Honky dory
-domain.set_quantity('friction', 100)     #Creep
+domain.set_quantity('friction', 100.0)     #Creep
 domain.set_quantity('elevation', island)
 domain.set_quantity('stage', 1)
@@ -94 +96 @@
 domain.check_integrity()
 
-
+domain.initialise_visualiser()
 #------------------------------------------------------------------------------
 # Evolve system through time

anuga_core/source/anuga/mesh_engine/test_triangmodule.py

@@ -211 +211 @@
 
     def testbad_hole(self):
-
+        
         holelist = [(9.0)]
         regionlist = []
-
-        points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-        pointattlist = [[],[],[],[]]
-        seglist = [(0,1),(1,3),(3,2),(2,0)]
+        
+        points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
+        pointattlist = [[],[],[],[]]
+        seglist = [(0,1),(1,3),(3,2),(2,0)]
         segattlist = [0,0,0,0]
         trilist = []

anuga_core/source/anuga/pmesh/mesh_interface.py

@@ -25 +25 @@
                              filename=None,
                              interior_regions=None,
+                             interior_holes=None,
                              poly_geo_reference=None,
                              mesh_geo_reference=None,
@@ -43 +44 @@
     Interior_regions is a list of tuples consisting of (polygon, resolution)
     for each region to be separately refined.
+    
+    Interior_holes is a list of ploygons for each hole.
 
     This function does not allow segments to share points - use underlying
@@ -71 +74 @@
               'filename': filename,
               'interior_regions': interior_regions,
+              'interior_holes': interior_holes,
               'poly_geo_reference': poly_geo_reference,
               'mesh_geo_reference': mesh_geo_reference,
@@ -106 +110 @@
                               filename=None,
                               interior_regions=None,
+                              interior_holes=None,
                               poly_geo_reference=None,
                               mesh_geo_reference=None,
@@ -146 +151 @@
                 raise PolygonError, msg
 
-
+    if interior_holes is not None:        
+        # Test that all the interior polygons are inside the bounding_poly
+        for interior_polygon, res in interior_holes:
+            indices = inside_polygon(interior_polygon, bounding_polygon,
+                                     closed = True, verbose = False)
+    
+            if len(indices) <> len(interior_polygon): 
+                msg = 'Interior polygon %s is outside bounding polygon: %s'\
+                      %(str(interior_polygon), str(bounding_polygon))
+                raise PolygonError, msg
 
     # Resolve geo referencing        
@@ -192 +206 @@
 
     # Do interior regions
+#    if interior_regions is not None:    
+#        for polygon, res in interior_regions:
+#            m.add_region_from_polygon(polygon,
+#                                      geo_reference=poly_geo_reference)
+#            # convert bounding poly to absolute values
+#            if poly_geo_reference is not None:
+#                polygon_absolute = \
+#                    poly_geo_reference.get_absolute(polygon)
+#            else:
+#                polygon_absolute = polygon
+#            inner_point = point_in_polygon(polygon_absolute)
+#            region = m.add_region(inner_point[0], inner_point[1])
+#            region.setMaxArea(res)
+            
+            
     if interior_regions is not None:    
         for polygon, res in interior_regions:
             m.add_region_from_polygon(polygon,
+                                      max_triangle_area = res,
                                       geo_reference=poly_geo_reference)
-            # convert bounding poly to absolute values
-            if poly_geo_reference is not None:
-                polygon_absolute = \
-                    poly_geo_reference.get_absolute(polygon)
-            else:
-                polygon_absolute = polygon
-            inner_point = point_in_polygon(polygon_absolute)
-            region = m.add_region(inner_point[0], inner_point[1])
-            region.setMaxArea(res)
+    
+            
+   # Do interior holes
+    if interior_holes is not None:    
+        for polygon, res in interior_holes:
+            m.add_hole_from_polygon(polygon,
+                                      geo_reference=poly_geo_reference)
+       
             
 

anuga_core/source/anuga/pmesh/test_mesh_interface.py

@@ -137 +137 @@
         interior_regions = [(inner1_polygon, 5),(inner2_polygon, 0.2)]
 
-
+        interior_holes = None
 
         # Clear cache first
@@ -147 +147 @@
                'filename': None,
                'interior_regions': interior_regions,
+               'interior_holes': interior_holes,
                'poly_geo_reference': geo_ref_poly,
                'mesh_geo_reference': mesh_geo,

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -524 +524 @@
 
     if (hc < minimum_allowed_height) {
-
+        
       //Old code: Set momentum to zero and ensure h is non negative
       //xmomc[k] = 0.0;
@@ -536 +536 @@
       //             has no Python equivalent
       if (hc <= 0.0) {
-        wc[k] = zc[k];
-        xmomc[k] = 0.0;
-        ymomc[k] = 0.0;
+        wc[k] = zc[k];
+                xmomc[k] = 0.0;
+                ymomc[k] = 0.0;
       } else {
         //Reduce excessive speeds derived from division by small hc
+        
         u = xmomc[k]/hc;
-        if (fabs(u) > maximum_allowed_speed) {
-          reduced_speed = maximum_allowed_speed * u/fabs(u);
-          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
-          //     u, reduced_speed);
-          xmomc[k] = reduced_speed * hc;
-        }
+                if (fabs(u) > maximum_allowed_speed) {
+                        reduced_speed = maximum_allowed_speed * u/fabs(u);
+                        //printf("Speed (u) has been reduced from %.3f to %.3f\n",
+                        //       u, reduced_speed);
+                        xmomc[k] = reduced_speed * hc;
+                }
 
         v = ymomc[k]/hc;
-        if (fabs(v) > maximum_allowed_speed) {
-          reduced_speed = maximum_allowed_speed * v/fabs(v);
-          //printf("Speed (v) has been reduced from %.3f to %.3f\n",
-          //     v, reduced_speed);
-          ymomc[k] = reduced_speed * hc;
-        }
+                if (fabs(v) > maximum_allowed_speed) {
+                        reduced_speed = maximum_allowed_speed * v/fabs(v);
+                        //printf("Speed (v) has been reduced from %.3f to %.3f\n",
+                        //       v, reduced_speed);
+                        ymomc[k] = reduced_speed * hc;
+                }
       }
     }