Changeset 6244

Timestamp:

Jan 30, 2009, 12:51:42 PM (16 years ago)

Author:

ole

Message:

Added better error messages and a new unit test that reveals
problem where fitting crashes by not being able to find
a triangle for a particular point within the mesh.

The test passes now, but if flag controlling mesh reuse
in quantity.py is switched on, it fails (on nautilus).

See ticket:314 for more info.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (1 diff)
• fit_interpolate/fit.py (modified) (6 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -874 +874 @@
             raise Exception, msg
 
-        # FIXME(Ole): I noticed a couple of examplse where this caused
+            
+        # FIXME(Ole): I noticed a couple of examples where this caused
         # a crash in fittng, so disabled it until I can investigate further
         # Sorry. 23 Jan 2009. Logged as ticket:314

anuga_core/source/anuga/fit_interpolate/fit.py

@@ -34 +34 @@
 from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
 from anuga.utilities.sparse import Sparse, Sparse_CSR
-from anuga.utilities.polygon import inside_polygon
+from anuga.utilities.polygon import inside_polygon, is_inside_polygon
 from anuga.fit_interpolate.search_functions import search_tree_of_vertices
 
@@ -115 +115 @@
             self._build_smoothing_matrix_D()
             
-        self.mesh_boundary_polygon = self.mesh.get_boundary_polygon()    
+        bd_poly = self.mesh.get_boundary_polygon()    
+        self.mesh_boundary_polygon = ensure_numeric(bd_poly)
             
     def _build_coefficient_matrix_B(self,
@@ -277 +278 @@
             # if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
             x = point_coordinates[i]
+            
             element_found, sigma0, sigma1, sigma2, k = \
                            search_tree_of_vertices(self.root, self.mesh, x)
@@ -300 +302 @@
             else:
                 # FIXME(Ole): This is the message referred to in ticket:314
+                
+                flag = is_inside_polygon(x,
+                                         self.mesh_boundary_polygon,
+                                         closed=True,
+                                         verbose=False) # Too much output if True                
+                msg = 'Point (%f, %f) is not inside mesh boundary' % tuple(x)
+                assert flag is True, msg                
+                
+                minx = min(self.mesh_boundary_polygon[:,0])
+                maxx = max(self.mesh_boundary_polygon[:,0])                
+                miny = min(self.mesh_boundary_polygon[:,1])
+                maxy = max(self.mesh_boundary_polygon[:,1])
                 msg = 'Could not find triangle for point %s. ' % str(x) 
-                msg += 'Mesh boundary is %s' % str(self.mesh_boundary_polygon)
+                msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
+                    % (minx, maxx, miny, maxy)
+                
+
                 raise Exception(msg)
             self.AtA = AtA
@@ -572 +589 @@
     """
 
-    # Duncan and Ole think that this isn't worth caching.
-    # Caching happens at the higher level anyway.
-    
-
     if mesh is None:
         # FIXME(DSG): Throw errors if triangles or vertex_coordinates
@@ -588 +601 @@
         mesh = Mesh(vertex_coordinates, triangles)
         mesh.check_integrity()
+    
     
     interp = Fit(mesh=mesh,

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -6526 +6526 @@
         os.remove(ptsfile)
 
+    def test_fitting_example_that_crashed(self):
+        """test_fitting_example_that_crashed
+        
+        This unit test has been derived from a real world example (the Towradgi '98 rainstorm simulation).
+        
+        It shows a condition where fitting as called from set_quantity crashes when ANUGA mesh
+        is reused. The test passes in the case where a new mesh is created.
+        
+        See ticket:314
+        """
+
+        verbose = False        
+
+        from anuga.shallow_water import Domain
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+        from anuga.geospatial_data.geospatial_data import Geospatial_data
+
+
+        #------------------------------------------------------------------------------
+        # Create domain
+        #------------------------------------------------------------------------------
+
+        W=303400
+        N=6195800
+        E=308640
+        S=6193120
+        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
+
+
+        offending_regions = []
+        # From culvert 8
+        offending_regions.append([[307611.43896231, 6193631.6894806],
+                                 [307600.11394969, 6193608.2855474],
+                                 [307597.41349586, 6193609.59227963],
+                                 [307608.73850848, 6193632.99621282]])
+        offending_regions.append([[307633.69143231, 6193620.9216536],
+                                 [307622.36641969, 6193597.5177204],
+                                 [307625.06687352, 6193596.21098818],
+                                 [307636.39188614, 6193619.61492137]])
+        # From culvert 9
+        offending_regions.append([[306326.69660524, 6194818.62900522],
+                                 [306324.67939476, 6194804.37099478],
+                                 [306323.75856492, 6194804.50127295],
+                                 [306325.7757754, 6194818.7592834]])
+        offending_regions.append([[306365.57160524, 6194813.12900522],
+                                 [306363.55439476, 6194798.87099478],
+                                 [306364.4752246, 6194798.7407166],
+                                 [306366.49243508, 6194812.99872705]])
+        # From culvert 10                         
+        offending_regions.append([[306955.071019428608, 6194465.704096679576],
+                                 [306951.616980571358, 6194457.295903320424],
+                                 [306950.044491164153, 6194457.941873183474],
+                                 [306953.498530021403, 6194466.350066542625]])
+        offending_regions.append([[307002.540019428649, 6194446.204096679576],
+                                 [306999.085980571399, 6194437.795903320424],
+                                 [307000.658469978604, 6194437.149933457375],
+                                 [307004.112508835853, 6194445.558126816526]])
+
+        interior_regions = []
+        for polygon in offending_regions:
+            interior_regions.append( [polygon, 100] ) 
+
+        meshname = 'offending_mesh.msh'
+        create_mesh_from_regions(bounding_polygon,
+                                 boundary_tags={'south': [0], 'east': [1], 'north': [2], 'west': [3]},
+                                 maximum_triangle_area=1000000,
+                                 interior_regions=interior_regions,
+                                 filename=meshname,
+                                 use_cache=False,
+                                 verbose=verbose)
+
+        domain = Domain(meshname, use_cache=False, verbose=verbose)
+
+
+        #------------------------------------------------------------------------------
+        # Fit data point to mesh
+        #------------------------------------------------------------------------------
+
+        points_file = 'offending_point.pts'
+
+        G=Geospatial_data(data_points=[[306953.344, 6194461.5]], # Offending point
+                          attributes=[1])
+        G.export_points_file(points_file)
+
+        domain.set_quantity('elevation', 
+                            filename=points_file,
+                            use_cache=False,
+                            verbose=verbose,
+                            alpha=0.01)
+
+