Changeset 5321


Ignore:
Timestamp:
May 14, 2008, 2:25:10 PM (16 years ago)
Author:
ole
Message:

Investigated performance profile of get_flow_through_cross_sections.
Most time appear to be spent in the interpolation routines.

Location:
anuga_core/source/anuga
Files:
6 edited

Legend:

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

    r5288 r5321  
    904904
    905905
    906     def _get_intersecting_segments(self, line):
     906    def _get_intersecting_segments(self, line,
     907                                   verbose=False):
    907908      """Find edges intersected by line
    908909
    909910      Input:
    910           line - list of two points forming a segmented line 
     911          line - list of two points forming a segmented line
     912          verbose
    911913      Output:
    912914          list of instances of class Triangle_intersection
     
    10471049
    10481050
    1049     def get_intersecting_segments(self, polyline):
     1051    def get_intersecting_segments(self, polyline,
     1052                                  verbose=False):
    10501053      """Find edges intersected by polyline
    10511054
    10521055      Input:
    10531056          polyline - list of points forming a segmented line
     1057          verbose
    10541058
    10551059      Output:
     
    10761080      triangle_intersections = []
    10771081      for i, point0 in enumerate(polyline[:-1]):
     1082
    10781083          point1 = polyline[i+1]
     1084          if verbose:
     1085              print 'Extracting mesh intersections from line:',
     1086              print '(%.2f, %.2f) - (%.2f, %.2f)' %(point0[0], point0[1],
     1087                                                    point1[0], point1[1])
     1088             
    10791089         
    10801090          line = [point0, point1]
  • anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py

    r4859 r5321  
    102102            if verbose: print 'FitInterpolate: Building mesh'       
    103103            self.mesh = Mesh(vertex_coordinates, triangles)
    104             self.mesh.check_integrity()
     104            #self.mesh.check_integrity() # Time consuming
    105105        else:
    106106            self.mesh = mesh
  • anuga_core/source/anuga/fit_interpolate/interpolate.py

    r5288 r5321  
    123123        # method is called even if interpolation points are unchanged.
    124124
    125 
    126         #print "point_coordinates interpolate.interpolate", point_coordinates
     125        #print "point_coordinates interpolate.interpolate", point_coordinates
     126        if verbose: print 'Build intepolation object'
    127127        if isinstance(point_coordinates, Geospatial_data):
    128128            point_coordinates = point_coordinates.get_data_points( \
     
    144144                msg = 'ERROR (interpolate.py): No point_coordinates inputted'
    145145                raise Exception(msg)
    146 
    147        
    148         if point_coordinates is not None:   
     146           
     147        if point_coordinates is not None:
    149148            self._point_coordinates = point_coordinates
    150149            if len(point_coordinates) < start_blocking_len or \
     
    202201        f = ensure_numeric(f, Float)       
    203202           
    204 
    205203        self._A = self._build_interpolation_matrix_A(point_coordinates,
    206204                                                     verbose=verbose)
     
    262260        """
    263261
    264         #print 'Building interpolation matrix'
     262        if verbose: print 'Building interpolation matrix'
    265263
    266264        # Convert point_coordinates to Numeric arrays, in case it was a list.
     
    670668                    result = interpol.interpolate(Q,
    671669                                                  point_coordinates=\
    672                                                   self.interpolation_points)
     670                                                  self.interpolation_points,
     671                                                  verbose=False) # Don't clutter
    673672
    674673                    #assert len(result), len(interpolation_points)
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5307 r5321  
    56555655   
    56565656    # Find all intersections and associated triangles.
    5657     segments = mesh.get_intersecting_segments(polyline)
     5657    segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
    56585658    #print
    56595659    #for x in segments:
     
    56685668        interpolation_points.append(midpoint)
    56695669
     5670
     5671    if verbose:
     5672        print 'Interpolating - ',       
     5673        print 'total number of interpolation points = %d'\
     5674              %len(interpolation_points)
    56705675
    56715676    I = Interpolation_function(time,
     
    56775682                               verbose=verbose)
    56785683
     5684    if verbose: print 'Computing hydrograph'
    56795685    # Compute hydrograph
    56805686    Q = []
  • anuga_core/source/anuga/utilities/polygon.py

    r5286 r5321  
    106106            if point_on_line([x3, y3], line0):
    107107                line1_ends_on_line0 = True                               
    108 
    109             #print line0_starts_on_line1
    110             #print line1_starts_on_line0           
    111             #print line1_ends_on_line0
    112             #print line0_ends_on_line1           
    113108
    114109            if not(line0_starts_on_line1 or line0_ends_on_line1\
     
    173168            # No intersection
    174169            return 0, None
     170
     171
     172def NEW_C_intersection(line0, line1):
     173    #FIXME(Ole): To write in C
     174    """Returns intersecting point between two line segments or None
     175    (if parallel or no intersection is found).
     176
     177    However, if parallel lines coincide partly (i.e. shara a common segment,
     178    the line segment where lines coincide is returned
     179   
     180
     181    Inputs:
     182        line0, line1: Each defined by two end points as in: [[x0, y0], [x1, y1]]
     183                      A line can also be a 2x2 numeric array with each row
     184                      corresponding to a point.
     185
     186
     187    Output:
     188        status, value
     189
     190        where status is interpreted as follows
     191       
     192        status == 0: no intersection with value set to None
     193        status == 1: One intersection point found and returned in value as [x,y]
     194        status == 2: Coinciding line segment found. Value taks the form [[x0,y0], [x1,y1]]
     195        status == 3: Lines would coincide but only if extended. Value set to None
     196        status == 4: Lines are parallel with a fixed distance apart. Value set to None.
     197   
     198    """
     199
     200
     201    line0 = ensure_numeric(line0, Float)
     202    line1 = ensure_numeric(line1, Float)   
     203
     204    status, value = _intersection(line0[0,0], line0[0,1],
     205                                  line0[1,0], line0[1,1],
     206                                  line1[0,0], line1[0,1],
     207                                  line1[1,0], line1[1,1])
     208
     209    return status, value
     210
    175211
    176212
     
    914950    from polygon_ext import _point_on_line
    915951    from polygon_ext import _separate_points_by_polygon
     952    #from polygon_ext import _intersection
    916953
    917954else:
  • anuga_core/source/anuga/utilities/polygon_ext.c

    r5225 r5321  
    5858}
    5959
     60
     61
     62/*
     63WORK IN PROGRESS TO OPTIMISE INTERSECTION
     64int __intersection(double x0, double y0,
     65                   double x1, double y1) {
     66
     67
     68    x0 = line0[0,0]; y0 = line0[0,1]
     69    x1 = line0[1,0]; y1 = line0[1,1]
     70
     71    x2 = line1[0,0]; y2 = line1[0,1]
     72    x3 = line1[1,0]; y3 = line1[1,1]
     73
     74    denom = (y3-y2)*(x1-x0) - (x3-x2)*(y1-y0)
     75    u0 = (x3-x2)*(y0-y2) - (y3-y2)*(x0-x2)
     76    u1 = (x2-x0)*(y1-y0) - (y2-y0)*(x1-x0)
     77       
     78    if allclose(denom, 0.0):
     79        # Lines are parallel - check if they coincide on a shared a segment
     80
     81        if allclose( [u0, u1], 0.0 ):
     82            # We now know that the lines if continued coincide
     83            # The remaining check will establish if the finite lines share a segment
     84
     85            line0_starts_on_line1 = line0_ends_on_line1 =\
     86            line1_starts_on_line0 = line1_ends_on_line0 = False
     87               
     88            if point_on_line([x0, y0], line1):
     89                line0_starts_on_line1 = True
     90
     91            if point_on_line([x1, y1], line1):
     92                line0_ends_on_line1 = True
     93 
     94            if point_on_line([x2, y2], line0):
     95                line1_starts_on_line0 = True
     96
     97            if point_on_line([x3, y3], line0):
     98                line1_ends_on_line0 = True                               
     99
     100            if not(line0_starts_on_line1 or line0_ends_on_line1\
     101               or line1_starts_on_line0 or line1_ends_on_line0):
     102                # Lines are parallel and would coincide if extended, but not as they are.
     103                return 3, None
     104
     105
     106            # One line fully included in the other. Use direction of included line
     107            if line0_starts_on_line1 and line0_ends_on_line1:
     108                # Shared segment is line0 fully included in line1
     109                segment = array([[x0, y0], [x1, y1]])               
     110
     111            if line1_starts_on_line0 and line1_ends_on_line0:
     112                # Shared segment is line1 fully included in line0
     113                segment = array([[x2, y2], [x3, y3]])
     114           
     115
     116            # Overlap with lines are oriented the same way
     117            if line0_starts_on_line1 and line1_ends_on_line0:
     118                # Shared segment from line0 start to line 1 end
     119                segment = array([[x0, y0], [x3, y3]])
     120
     121            if line1_starts_on_line0 and line0_ends_on_line1:
     122                # Shared segment from line1 start to line 0 end
     123                segment = array([[x2, y2], [x1, y1]])                               
     124
     125
     126            # Overlap in opposite directions - use direction of line0
     127            if line0_starts_on_line1 and line1_starts_on_line0:
     128                # Shared segment from line0 start to line 1 end
     129                segment = array([[x0, y0], [x2, y2]])
     130
     131            if line0_ends_on_line1 and line1_ends_on_line0:
     132                # Shared segment from line0 start to line 1 end
     133                segment = array([[x3, y3], [x1, y1]])               
     134
     135               
     136            return 2, segment
     137        else:
     138            # Lines are parallel but they do not coincide
     139            return 4, None #FIXME (Ole): Add distance here instead of None
     140           
     141    else:
     142        # Lines are not parallel or coinciding
     143        u0 = u0/denom
     144        u1 = u1/denom       
     145
     146        x = x0 + u0*(x1-x0)
     147        y = y0 + u0*(y1-y0)
     148
     149        # Sanity check - can be removed to speed up if needed
     150        assert allclose(x, x2 + u1*(x3-x2))
     151        assert allclose(y, y2 + u1*(y3-y2))       
     152
     153        # Check if point found lies within given line segments
     154        if 0.0 <= u0 <= 1.0 and 0.0 <= u1 <= 1.0:
     155            # We have intersection
     156
     157            return 1, array([x, y])
     158        else:
     159            # No intersection
     160            return 0, None
     161
     162
     163}
     164*/
    60165
    61166
     
    176281
    177282
     283/*
     284PyObject *_intersection(PyObject *self, PyObject *args) {
     285  //
     286  // intersection(x0, y0, x1, y1)
     287  //
     288
     289  double x, y, x0, y0, x1, y1;
     290  int res;
     291  PyObject *result;
     292
     293  // Convert Python arguments to C
     294  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) {
     295    PyErr_SetString(PyExc_RuntimeError,
     296                    "point_on_line could not parse input");   
     297    return NULL;
     298  }
     299
     300
     301  // Call underlying routine
     302  res = __intersection(x, y, x0, y0, x1, y1);
     303
     304  // Return values a and b
     305  result = Py_BuildValue("i", res);
     306  return result;
     307}
     308*/
     309
    178310
    179311PyObject *_separate_points_by_polygon(PyObject *self, PyObject *args) {
     
    266398
    267399  {"_point_on_line", _point_on_line, METH_VARARGS, "Print out"},
     400  //{"_intersection", _intersection, METH_VARARGS, "Print out"}, 
    268401  {"_separate_points_by_polygon", _separate_points_by_polygon,
    269402                                 METH_VARARGS, "Print out"},
Note: See TracChangeset for help on using the changeset viewer.