Changeset 6236


Ignore:
Timestamp:
Jan 29, 2009, 2:37:50 AM (10 years ago)
Author:
ole
Message:

Better error message in Fit and cleanup.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/fit.py

    r6232 r6236  
    3434from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate
    3535from anuga.utilities.sparse import Sparse, Sparse_CSR
    36 from anuga.utilities.polygon import in_and_outside_polygon
     36from anuga.utilities.polygon import inside_polygon
    3737from anuga.fit_interpolate.search_functions import search_tree_of_vertices
    3838
     
    4646
    4747import Numeric as num
    48 
    49 
    50 #DEFAULT_ALPHA = 0.001
    5148
    5249
     
    137134            self.B = self.AtA
    138135
    139         #Convert self.B matrix to CSR format for faster matrix vector
     136        # Convert self.B matrix to CSR format for faster matrix vector
    140137        self.B = Sparse_CSR(self.B)
    141138
     
    166163        """
    167164       
    168         #FIXME: algorithm might be optimised by computing local 9x9
    169         #"element stiffness matrices:
     165        # FIXME: algorithm might be optimised by computing local 9x9
     166        # "element stiffness matrices:
    170167
    171168        m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
     
    173170        self.D = Sparse(m,m)
    174171
    175         #For each triangle compute contributions to D = D1+D2
     172        # For each triangle compute contributions to D = D1+D2
    176173        for i in range(len(self.mesh)):
    177174
    178             #Get area
     175            # Get area
    179176            area = self.mesh.areas[i]
    180177
    181             #Get global vertex indices
     178            # Get global vertex indices
    182179            v0 = self.mesh.triangles[i,0]
    183180            v1 = self.mesh.triangles[i,1]
    184181            v2 = self.mesh.triangles[i,2]
    185182
    186             #Get the three vertex_points
     183            # Get the three vertex_points
    187184            xi0 = self.mesh.get_vertex_coordinate(i, 0)
    188185            xi1 = self.mesh.get_vertex_coordinate(i, 1)
    189186            xi2 = self.mesh.get_vertex_coordinate(i, 2)
    190187
    191             #Compute gradients for each vertex
     188            # Compute gradients for each vertex
    192189            a0, b0 = gradient(xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1],
    193190                              1, 0, 0)
     
    199196                              0, 0, 1)
    200197
    201             #Compute diagonal contributions
     198            # Compute diagonal contributions
    202199            self.D[v0,v0] += (a0*a0 + b0*b0)*area
    203200            self.D[v1,v1] += (a1*a1 + b1*b1)*area
    204201            self.D[v2,v2] += (a2*a2 + b2*b2)*area
    205202
    206             #Compute contributions for basis functions sharing edges
     203            # Compute contributions for basis functions sharing edges
    207204            e01 = (a0*a1 + b0*b1)*area
    208205            self.D[v0,v1] += e01
     
    246243        The number of attributes of the data points does not change
    247244        """
    248         #Build n x m interpolation matrix
    249 
     245       
     246        # Build n x m interpolation matrix
    250247        if self.AtA == None:
    251248            # AtA and Atz need to be initialised.
     
    265262        self.point_count += point_coordinates.shape[0]
    266263
    267         #if verbose: print 'Getting indices inside mesh boundary'
    268 
    269         inside_poly_indices, outside_poly_indices  = \
    270                      in_and_outside_polygon(point_coordinates,
    271                                             self.mesh_boundary_polygon,
    272                                             closed = True,
    273                                             verbose = False) # There's too much output if True
    274         #print "self.inside_poly_indices",self.inside_poly_indices
    275         #print "self.outside_poly_indices",self.outside_poly_indices
    276 
    277        
    278         n = len(inside_poly_indices)
    279         #if verbose: print 'Building fitting matrix from %d points' %n       
    280 
    281         #Compute matrix elements for points inside the mesh
    282         triangles = self.mesh.triangles #Did this for speed, did ~nothing
    283         for d, i in enumerate(inside_poly_indices):
     264
     265        inside_indices = inside_polygon(point_coordinates,
     266                                        self.mesh_boundary_polygon,
     267                                        closed=True,
     268                                        verbose=False) # Too much output if True
     269
     270       
     271        n = len(inside_indices)
     272
     273        # Compute matrix elements for points inside the mesh
     274        triangles = self.mesh.triangles # Shorthand
     275        for d, i in enumerate(inside_indices):
    284276            # For each data_coordinate point
    285277            # if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n)
     
    307299                        AtA[j,k] += sigmas[j]*sigmas[k]
    308300            else:
    309                 msg = 'Could not find triangle for point', x
     301                msg = 'Could not find triangle for point %s. ' % str(x)
     302                msg += 'Mesh boundary is %s' % str(self.mesh_boundary_polygon)
    310303                raise Exception(msg)
    311304            self.AtA = AtA
     305
    312306       
    313307    def fit(self, point_coordinates_or_filename=None, z=None,
     
    329323         
    330324        """
    331         # use blocking to load in the point info
     325       
     326        # Use blocking to load in the point info
    332327        if type(point_coordinates_or_filename) == types.StringType:
    333328            msg = "Don't set a point origin when reading from a file"
     
    359354
    360355                points = geo_block.get_data_points(absolute=True)
    361                 #print "fit points", points
    362356                z = geo_block.get_attributes(attribute_name=attribute_name)
    363357                self.build_fit_subset(points, z, verbose=verbose)
     
    373367            assert self.Atz <> None
    374368           
    375             #FIXME (DSG) - do  a message
     369            # FIXME (DSG) - do  a message
    376370        else:
    377371            point_coordinates = ensure_absolute(point_coordinates,
    378372                                                geo_reference=point_origin)
    379             #if isinstance(point_coordinates,Geospatial_data) and z is None:
     373            # if isinstance(point_coordinates,Geospatial_data) and z is None:
    380374            # z will come from the geo-ref
    381375            self.build_fit_subset(point_coordinates, z, verbose)
    382376
    383         #Check sanity
     377        # Check sanity
    384378        m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex)
    385379        n = self.point_count
     
    457451                alpha=DEFAULT_ALPHA,
    458452                verbose=False,
    459                 acceptable_overshoot=1.01, # FIXME: Move to config - this value is assumed in caching test
     453                acceptable_overshoot=1.01,
     454                # FIXME: Move to config - this value is assumed in caching test
     455                # FIXME(Ole): Just realised that this was never implemented (29 Jan 2009). I suggest removing it altogether.
    460456                mesh_origin=None,
    461457                data_origin=None,
Note: See TracChangeset for help on using the changeset viewer.