Changeset 6244


Ignore:
Timestamp:
Jan 30, 2009, 12:51:42 PM (10 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:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r6237 r6244  
    874874            raise Exception, msg
    875875
    876         # FIXME(Ole): I noticed a couple of examplse where this caused
     876           
     877        # FIXME(Ole): I noticed a couple of examples where this caused
    877878        # a crash in fittng, so disabled it until I can investigate further
    878879        # Sorry. 23 Jan 2009. Logged as ticket:314
  • anuga_core/source/anuga/fit_interpolate/fit.py

    r6237 r6244  
    3434from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
    3535from anuga.utilities.sparse import Sparse, Sparse_CSR
    36 from anuga.utilities.polygon import inside_polygon
     36from anuga.utilities.polygon import inside_polygon, is_inside_polygon
    3737from anuga.fit_interpolate.search_functions import search_tree_of_vertices
    3838
     
    115115            self._build_smoothing_matrix_D()
    116116           
    117         self.mesh_boundary_polygon = self.mesh.get_boundary_polygon()   
     117        bd_poly = self.mesh.get_boundary_polygon()   
     118        self.mesh_boundary_polygon = ensure_numeric(bd_poly)
    118119           
    119120    def _build_coefficient_matrix_B(self,
     
    277278            # if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
    278279            x = point_coordinates[i]
     280           
    279281            element_found, sigma0, sigma1, sigma2, k = \
    280282                           search_tree_of_vertices(self.root, self.mesh, x)
     
    300302            else:
    301303                # FIXME(Ole): This is the message referred to in ticket:314
     304               
     305                flag = is_inside_polygon(x,
     306                                         self.mesh_boundary_polygon,
     307                                         closed=True,
     308                                         verbose=False) # Too much output if True               
     309                msg = 'Point (%f, %f) is not inside mesh boundary' % tuple(x)
     310                assert flag is True, msg               
     311               
     312                minx = min(self.mesh_boundary_polygon[:,0])
     313                maxx = max(self.mesh_boundary_polygon[:,0])               
     314                miny = min(self.mesh_boundary_polygon[:,1])
     315                maxy = max(self.mesh_boundary_polygon[:,1])
    302316                msg = 'Could not find triangle for point %s. ' % str(x)
    303                 msg += 'Mesh boundary is %s' % str(self.mesh_boundary_polygon)
     317                msg += 'Mesh boundary extent is (%.f, %.f), (%.f, %.f)'\
     318                    % (minx, maxx, miny, maxy)
     319               
     320
    304321                raise Exception(msg)
    305322            self.AtA = AtA
     
    572589    """
    573590
    574     # Duncan and Ole think that this isn't worth caching.
    575     # Caching happens at the higher level anyway.
    576    
    577 
    578591    if mesh is None:
    579592        # FIXME(DSG): Throw errors if triangles or vertex_coordinates
     
    588601        mesh = Mesh(vertex_coordinates, triangles)
    589602        mesh.check_integrity()
     603   
    590604   
    591605    interp = Fit(mesh=mesh,
  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r6191 r6244  
    65266526        os.remove(ptsfile)
    65276527
     6528    def test_fitting_example_that_crashed(self):
     6529        """test_fitting_example_that_crashed
     6530       
     6531        This unit test has been derived from a real world example (the Towradgi '98 rainstorm simulation).
     6532       
     6533        It shows a condition where fitting as called from set_quantity crashes when ANUGA mesh
     6534        is reused. The test passes in the case where a new mesh is created.
     6535       
     6536        See ticket:314
     6537        """
     6538
     6539        verbose = False       
     6540
     6541        from anuga.shallow_water import Domain
     6542        from anuga.pmesh.mesh_interface import create_mesh_from_regions
     6543        from anuga.geospatial_data.geospatial_data import Geospatial_data
     6544
     6545
     6546        #------------------------------------------------------------------------------
     6547        # Create domain
     6548        #------------------------------------------------------------------------------
     6549
     6550        W=303400
     6551        N=6195800
     6552        E=308640
     6553        S=6193120
     6554        bounding_polygon = [[W, S], [E, S], [E, N], [W, N]]
     6555
     6556
     6557        offending_regions = []
     6558        # From culvert 8
     6559        offending_regions.append([[307611.43896231, 6193631.6894806],
     6560                                 [307600.11394969, 6193608.2855474],
     6561                                 [307597.41349586, 6193609.59227963],
     6562                                 [307608.73850848, 6193632.99621282]])
     6563        offending_regions.append([[307633.69143231, 6193620.9216536],
     6564                                 [307622.36641969, 6193597.5177204],
     6565                                 [307625.06687352, 6193596.21098818],
     6566                                 [307636.39188614, 6193619.61492137]])
     6567        # From culvert 9
     6568        offending_regions.append([[306326.69660524, 6194818.62900522],
     6569                                 [306324.67939476, 6194804.37099478],
     6570                                 [306323.75856492, 6194804.50127295],
     6571                                 [306325.7757754, 6194818.7592834]])
     6572        offending_regions.append([[306365.57160524, 6194813.12900522],
     6573                                 [306363.55439476, 6194798.87099478],
     6574                                 [306364.4752246, 6194798.7407166],
     6575                                 [306366.49243508, 6194812.99872705]])
     6576        # From culvert 10                         
     6577        offending_regions.append([[306955.071019428608, 6194465.704096679576],
     6578                                 [306951.616980571358, 6194457.295903320424],
     6579                                 [306950.044491164153, 6194457.941873183474],
     6580                                 [306953.498530021403, 6194466.350066542625]])
     6581        offending_regions.append([[307002.540019428649, 6194446.204096679576],
     6582                                 [306999.085980571399, 6194437.795903320424],
     6583                                 [307000.658469978604, 6194437.149933457375],
     6584                                 [307004.112508835853, 6194445.558126816526]])
     6585
     6586        interior_regions = []
     6587        for polygon in offending_regions:
     6588            interior_regions.append( [polygon, 100] )
     6589
     6590        meshname = 'offending_mesh.msh'
     6591        create_mesh_from_regions(bounding_polygon,
     6592                                 boundary_tags={'south': [0], 'east': [1], 'north': [2], 'west': [3]},
     6593                                 maximum_triangle_area=1000000,
     6594                                 interior_regions=interior_regions,
     6595                                 filename=meshname,
     6596                                 use_cache=False,
     6597                                 verbose=verbose)
     6598
     6599        domain = Domain(meshname, use_cache=False, verbose=verbose)
     6600
     6601
     6602        #------------------------------------------------------------------------------
     6603        # Fit data point to mesh
     6604        #------------------------------------------------------------------------------
     6605
     6606        points_file = 'offending_point.pts'
     6607
     6608        G=Geospatial_data(data_points=[[306953.344, 6194461.5]], # Offending point
     6609                          attributes=[1])
     6610        G.export_points_file(points_file)
     6611
     6612        domain.set_quantity('elevation',
     6613                            filename=points_file,
     6614                            use_cache=False,
     6615                            verbose=verbose,
     6616                            alpha=0.01)
     6617
     6618       
    65286619
    65296620       
Note: See TracChangeset for help on using the changeset viewer.