Changeset 1020


Ignore:
Timestamp:
Mar 7, 2005, 9:38:36 AM (20 years ago)
Author:
steve
Message:
 
Location:
inundation/ga/storm_surge/pyvolution
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/domain.py

    r991 r1020  
    2020        Mesh.__init__(self, coordinates, vertices, boundary, tagged_elements)
    2121
    22         from Numeric import zeros, Float 
     22        from Numeric import zeros, Float
    2323        from quantity import Quantity, Conserved_quantity
    2424
    … …  
    2626        #the conservation equations
    2727        #(Must be a subset of quantities)
    28         if conserved_quantities is None:           
     28        if conserved_quantities is None:
    2929            self.conserved_quantities = []
    3030        else:
    31             self.conserved_quantities = conserved_quantities           
     31            self.conserved_quantities = conserved_quantities
    3232
    3333        if other_quantities is None:
    3434            self.other_quantities = []
    3535        else:
    36             self.other_quantities = other_quantities           
    37            
     36            self.other_quantities = other_quantities
     37
    3838
    3939        #Build dictionary of Quantity instances keyed by quantity names
    … …  
    5353        from config import max_smallsteps, beta_w, beta_h, epsilon, CFL
    5454        self.beta_w = beta_w
    55         self.beta_h = beta_h       
     55        self.beta_h = beta_h
    5656        self.epsilon = epsilon
    5757
    … …  
    6161        self.order = self.default_order
    6262        self.smallsteps = 0
    63         self.max_smallsteps = max_smallsteps       
     63        self.max_smallsteps = max_smallsteps
    6464        self.number_of_steps = 0
    65         self.number_of_first_order_steps = 0       
     65        self.number_of_first_order_steps = 0
    6666        self.CFL = CFL
    67        
    68         #Model time   
     67
     68        #Model time
    6969        self.time = 0.0
    7070        self.finaltime = None
    … …  
    7474        self.zone = zone
    7575        self.xllcorner = xllcorner
    76         self.yllcorner = yllcorner       
    77        
    78        
    79         #Checkpointing and storage             
     76        self.yllcorner = yllcorner
     77
     78
     79        #Checkpointing and storage
    8080        from config import default_datadir
    81         self.datadir = default_datadir 
     81        self.datadir = default_datadir
    8282        self.filename = 'domain'
    8383        self.checkpoint = False
    84        
    85 
    86     #Public interface to Domain       
     84
     85
     86    #Public interface to Domain
    8787    def get_conserved_quantities(self, vol_id, vertex=None, edge=None):
    8888        """Get conserved quantities at volume vol_id
    … …  
    9393        If both are specified an exeception is raised
    9494
    95         Return value: Vector of length == number_of_conserved quantities 
    96        
     95        Return value: Vector of length == number_of_conserved quantities
     96
    9797        """
    9898
    9999        from Numeric import zeros, Float
    100        
     100
    101101        if not (vertex is None or edge is None):
    102102            msg = 'Values for both vertex and edge was specified.'
    103             msg += 'Only one (or none) is allowed.' 
     103            msg += 'Only one (or none) is allowed.'
    104104            raise msg
    105105
    106106        q = zeros( len(self.conserved_quantities), Float)
    107        
     107
    108108        for i, name in enumerate(self.conserved_quantities):
    109109            Q = self.quantities[name]
    … …  
    113113                q[i] = Q.edge_values[vol_id, edge]
    114114            else:
    115                 q[i] = Q.centroid_values[vol_id]               
    116 
    117         return q 
     115                q[i] = Q.centroid_values[vol_id]
     116
     117        return q
    118118
    119119
    … …  
    121121        """Set values for named quantities.
    122122        The index is the quantity
    123        
     123
    124124        name: Name of quantity
    125125        X: Compatible list, Numeric array, const or function (see below)
    126        
     126
    127127        The values will be stored in elements following their
    128128        internal ordering.
    … …  
    131131        for key in quantity_dict.keys():
    132132            self.set_quantity(key, quantity_dict[key], location='vertices')
    133        
     133
    134134
    135135    def set_quantity(self, name, X, location='vertices', indexes = None):
    136136        """Set values for named quantity
    137        
     137
    138138        name: Name of quantity
    139139        X: Compatible list, Numeric array, const or function (see below)
    … …  
    144144        be a list of a Numerical array of length N, N being the number
    145145        of elements. Otherwise it must be of dimension Nx3.
    146        
     146
    147147        Indexes is the set of element ids that the operation applies to
    148        
     148
    149149        The values will be stored in elements following their
    150150        internal ordering.
    151        
     151
    152152        #FIXME (Ole): I suggest the following interface
    153153        set_quantity(name, X, location, region)
    154154        where
    155155            name: Name of quantity
    156             X: 
    157               -Compatible list, 
    158               -Numeric array, 
     156            X:
     157              -Compatible list,
     158              -Numeric array,
    159159              -const or function (see below)
    160160              -another quantity Q or an expression of the form
    … …  
    168168                    - indices (refers to specific triangles)
    169169                    - polygon (identifies region)
    170                    
     170
    171171            This should work for all values of X
    172            
     172
    173173
    174174
    … …  
    176176
    177177        #from quantity import Quantity, Conserved_quantity
    178        
     178
    179179        #Create appropriate quantity object
    180180        ##if name in self.conserved_quantities:
    181181        ##    self.quantities[name] = Conserved_quantity(self)
    182182        ##else:
    183         ##    self.quantities[name] = Quantity(self)           
     183        ##    self.quantities[name] = Quantity(self)
    184184
    185185        #Set value
    186186        self.quantities[name].set_values(X, location, indexes = indexes)
    187187
    188          
     188
    189189    def get_quantity(self, name, location='vertices', indexes = None):
    190190        """Get values for named quantity
    191        
     191
    192192        name: Name of quantity
    193        
     193
    194194        In case of location == 'centroid' the dimension values must
    195195        be a list of a Numerical array of length N, N being the number
    … …  
    197197
    198198        Indexes is the set of element ids that the operation applies to.
    199        
     199
    200200        The values will be stored in elements following their
    201201        internal ordering.
    … …  
    203203
    204204        return self.quantities[name].get_values( location, indexes = indexes)
    205    
     205
    206206
    207207    def set_boundary(self, boundary_map):
    … …  
    215215        in the list self.boundary_objects.
    216216        More entries may point to the same boundary object
    217        
     217
    218218        Schematically the mapping is from two dictionaries to one list
    219219        where the index is used as pointer to the boundary_values arrays
    … …  
    224224        ----------------------------------------------
    225225        self.boundary_objects:  ((vol_id, edge_id), boundary_object)
    226        
     226
    227227
    228228        Pre-condition:
    … …  
    236236        However, if a tag is not used to the domain, no error is thrown.
    237237        FIXME: This would lead to implementation of a
    238         default boundary condition       
     238        default boundary condition
    239239
    240240        Note: If a segment is listed in the boundary dictionary and if it is
    … …  
    242242        even if there is a neighbouring triangle.
    243243        This would be the case for internal boundaries
    244        
     244
    245245        Boundary objects that are None will be skipped.
    246246
    … …  
    248248        object is changed into None, the neighbour structure will not be
    249249        restored!!!
    250        
     250
    251251
    252252        """
    … …  
    269269                    pass
    270270                    #FIXME: Check and perhaps fix neighbour structure
    271                    
     271
    272272
    273273            else:
    … …  
    279279                raise msg
    280280
    281                
     281
    282282    def set_region(self, functions):
    283283        # The order of functions in the list is used.
    … …  
    290290                #Do we need to do this sort of thing?
    291291                #self = function(tag, self.tagged_elements[tag], self)
    292            
    293     #MISC       
     292
     293    #MISC
    294294    def check_integrity(self):
    295295        Mesh.check_integrity(self)
    … …  
    299299            assert quantity in self.quantities, msg
    300300
    301         ##assert hasattr(self, 'boundary_objects')   
    302    
     301        ##assert hasattr(self, 'boundary_objects')
     302
    303303    def write_time(self):
    304304        if self.min_timestep == self.max_timestep:
    … …  
    309309            print 'Time = %.4f, steps=%d (%d)'\
    310310                  %(self.time, self.number_of_steps,
    311                     self.number_of_first_order_steps)           
     311                    self.number_of_first_order_steps)
    312312        else:
    313313            print 'Time = %.4f, delta t in [%.8f, %.8f], steps=%d (%d)'\
    … …  
    318318
    319319    def get_name(self):
    320         return self.filename   
    321    
     320        return self.filename
     321
    322322    def set_name(self, name):
    323         self.filename = name   
     323        self.filename = name
    324324
    325325    def get_datadir(self):
    326         return self.datadir     
    327                        
     326        return self.datadir
     327
    328328    def set_datadir(self, name):
    329         self.datadir = name     
    330 
    331        
     329        self.datadir = name
     330
     331
    332332
    333333    #def set_defaults(self):
    … …  
    337337    #
    338338    #    for name in self.conserved_quantities + self.other_quantities:
    339     #        self.set_quantity(name, 0.0)       
    340                
     339    #        self.set_quantity(name, 0.0)
     340
    341341
    342342    ###########################
    … …  
    346346        """Evolve model from time=0.0 to finaltime yielding results
    347347        every yieldstep.
    348        
     348
    349349        Internally, smaller timesteps may be taken.
    350350
    351351        Evolve is implemented as a generator and is to be called as such, e.g.
    352        
     352
    353353        for t in domain.evolve(timestep, yieldstep, finaltime):
    354354            <Do something with domain and t>
    355        
     355
    356356        """
    357357
    … …  
    360360        #FIXME: Maybe lump into a larger check prior to evolving
    361361        msg = 'Boundary tags must be bound to boundary objects before evolving system, '
    362         msg += 'e.g. using the method set_boundary.\n'       
     362        msg += 'e.g. using the method set_boundary.\n'
    363363        msg += 'This system has the boundary tags %s ' %self.get_boundary_tags()
    364364        assert hasattr(self, 'boundary_objects'), msg
    365        
     365
    366366        ##self.set_defaults()
    367367
    … …  
    370370
    371371        self.order = self.default_order
    372        
    373            
     372
     373
    374374        self.yieldtime = 0.0 #Time between 'yields'
    375375
    … …  
    379379        self.finaltime = finaltime
    380380        self.number_of_steps = 0
    381         self.number_of_first_order_steps = 0               
     381        self.number_of_first_order_steps = 0
    382382
    383383        #Initial update of vertex and edge values
    384384        self.distribute_to_vertices_and_edges()
    385        
    386        
     385
     386
    387387        #Or maybe restore from latest checkpoint
    388         if self.checkpoint is True:
     388        if self.checkpoint is True:
    389389            self.goto_latest_checkpoint()
    390390
    391            
     391
    392392        yield(self.time)  #Yield initial values
    393        
     393
    394394        while True:
    395395            #Update boundary values
    396             self.update_boundary()           
    397 
    398             #Compute fluxes across each element edge
     396            self.update_boundary()
     397
     398            #Compute fluxes across each element edge
    399399            self.compute_fluxes()
    400400
    … …  
    403403
    404404            #Update conserved quantities
    405             self.update_conserved_quantities()           
     405            self.update_conserved_quantities()
    406406
    407407            #Update vertex and edge values
    408408            self.distribute_to_vertices_and_edges()
    409            
    410             #Update time   
     409
     410            #Update time
    411411            self.time += self.timestep
    412412            self.yieldtime += self.timestep
    413413            self.number_of_steps += 1
    414414            if self.order == 1:
    415                 self.number_of_first_order_steps += 1       
     415                self.number_of_first_order_steps += 1
    416416
    417417            #Yield results
    … …  
    421421                break
    422422
    423                
    424             if abs(self.yieldtime - yieldstep) < epsilon: 
     423
     424            if abs(self.yieldtime - yieldstep) < epsilon:
    425425                # Yield (intermediate) time and allow inspection of domain
    426426
    … …  
    428428                    self.store_checkpoint()
    429429                    self.delete_old_checkpoints()
    430                    
    431                 #Pass control on to outer loop for more specific actions   
    432                 yield(self.time) 
     430
     431                #Pass control on to outer loop for more specific actions
     432                yield(self.time)
    433433
    434434                # Reinitialise
    435                 self.yieldtime = 0.0               
     435                self.yieldtime = 0.0
    436436                self.min_timestep = max_timestep
    437437                self.max_timestep = min_timestep
    438438                self.number_of_steps = 0
    439                 self.number_of_first_order_steps = 0                       
    440 
    441            
     439                self.number_of_first_order_steps = 0
     440
     441
    442442    def evolve_to_end(self, finaltime = 1.0):
    443443        """Iterate evolve all the way to the end
    444444        """
    445445
    446         for _ in self.evolve(yieldstep=None, finaltime=finaltime):           
     446        for _ in self.evolve(yieldstep=None, finaltime=finaltime):
    447447            pass
    448        
    449 
    450    
     448
     449
     450
    451451    def update_boundary(self):
    452452        """Go through list of boundary objects and update boundary values
    … …  
    471471
    472472        from config import min_timestep
    473    
     473
    474474        timestep = self.timestep
    475    
    476         #Record maximal and minimal values of timestep for reporting   
     475
     476        #Record maximal and minimal values of timestep for reporting
    477477        self.max_timestep = max(timestep, self.max_timestep)
    478478        self.min_timestep = min(timestep, self.min_timestep)
    479    
     479
    480480        #Protect against degenerate time steps
    481481        if timestep < min_timestep:
    482482
    483483            #Number of consecutive small steps taken b4 taking action
    484             self.smallsteps += 1 
     484            self.smallsteps += 1
    485485
    486486            if self.smallsteps > self.max_smallsteps:
    487487                self.smallsteps = 0 #Reset
    488                
     488
    489489                if self.order == 1:
    490490                    msg = 'Too small timestep %.16f reached ' %timestep
    … …  
    501501            if self.order == 1 and self.default_order == 2:
    502502                self.order = 2
    503                    
     503
    504504
    505505        #Ensure that final time is not exceeded
    … …  
    515515
    516516
    517     def compute_forcing_terms(self):         
     517    def compute_forcing_terms(self):
    518518        """If there are any forcing functions driving the system
    519519        they should be defined in Domain subclass and appended to
    … …  
    524524            f(self)
    525525
    526        
     526
    527527
    528528    def update_conserved_quantities(self):
    … …  
    532532
    533533        from Numeric import ones, sum, equal, Float
    534            
     534
    535535        N = self.number_of_elements
    536536        d = len(self.conserved_quantities)
    537        
     537
    538538        timestep = self.timestep
    539539
    540         #Compute forcing terms   
    541         self.compute_forcing_terms()   
     540        #Compute forcing terms
     541        self.compute_forcing_terms()
    542542
    543543        #Update conserved_quantities
    … …  
    547547
    548548            #Clean up
    549             #Note that Q.explicit_update is reset by compute_fluxes           
     549            #Note that Q.explicit_update is reset by compute_fluxes
    550550            Q.semi_implicit_update[:] = 0.0
    551551
    552            
     552
    553553
    554554    def distribute_to_vertices_and_edges(self):
    … …  
    567567            elif self.order == 2:
    568568                Q.extrapolate_second_order()
    569                 Q.limit()                               
     569                Q.limit()
    570570            else:
    571571                raise 'Unknown order'
    572             Q.interpolate_from_vertices_to_edges()                           
    573 
    574 
    575                
     572            Q.interpolate_from_vertices_to_edges()
     573
     574
     575
    576576##############################################
    577577#Initialise module
    … …  
    587587        print msg
    588588    else:
    589         psyco.bind(Domain.update_boundary) 
     589        psyco.bind(Domain.update_boundary)
    590590        #psyco.bind(Domain.update_timestep)     #Not worth it
    591591        psyco.bind(Domain.update_conserved_quantities)
    592592        psyco.bind(Domain.distribute_to_vertices_and_edges)
    593        
     593
    594594
    595595if __name__ == "__main__":

  • inundation/ga/storm_surge/pyvolution/util_ext.c

    r999 r1020  
    22//
    33// To compile (Python2.3):
    4 //  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O 
    5 //  gcc -shared util_ext.o  -o util_ext.so     
     4//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O
     5//  gcc -shared util_ext.o  -o util_ext.so
    66//
    77// See the module util.py
    … …  
    99//
    1010// Ole Nielsen, GA 2004
    11        
     11
    1212#include "Python.h"
    1313#include "Numeric/arrayobject.h"
    … …  
    1919
    2020
    21 int _point_on_line(double x, double y, 
    22                    double x0, double y0, 
     21int _point_on_line(double x, double y,
     22                   double x0, double y0,
    2323                   double x1, double y1) {
    24   /*Determine whether a point is on a line segment 
    25    
     24  /*Determine whether a point is on a line segment
     25
    2626    Input: x, y, x0, x0, x1, y1: where
    2727        point is given by x, y
    2828        line is given by (x0, y0) and (x1, y1)
    29        
     29
    3030  */
    31    
     31
    3232  double a0, a1, a_normal0, a_normal1, b0, b1, len_a, len_b;
    33  
     33
    3434  a0 = x - x0;
    35   a1 = y - y0; 
    36  
     35  a1 = y - y0;
     36
    3737  a_normal0 = a1;
    3838  a_normal1 = -a0;
    39                        
     39
    4040  b0 = x1 - x0;
    4141  b1 = y1 - y0;
    42    
    43   if ( a_normal0*b0 + a_normal1*b1 == 0 ) {   
     42
     43  if ( a_normal0*b0 + a_normal1*b1 == 0 ) {
    4444    //Point is somewhere on the infinite extension of the line
    4545
    46     len_a = sqrt(a0*a0 + a1*a1);               
     46    len_a = sqrt(a0*a0 + a1*a1);
    4747    len_b = sqrt(b0*b0 + b1*b1);
    48                                        
     48
    4949    if (a0*b0 + a1*b1 >= 0 && len_a <= len_b) {
    5050      return 1;
    … …  
    6363                    double* points,
    6464                    double* polygon,
    65                     int* indices,    // M-Array for storage indices 
    66                     int closed, 
     65                    int* indices,    // M-Array for storage indices
     66                    int closed,
    6767                    int verbose) {
    68  
     68
    6969  double minpx, maxpx, minpy, maxpy, x, y, px_i, py_i, px_j, py_j;
    7070  int i, j, k, count, inside;
    71  
    72   //Find min and max of poly used for optimisation when points 
    73   //are far away from polygon 
    74  
    75  
     71
     72  //Find min and max of poly used for optimisation when points
     73  //are far away from polygon
     74
     75
    7676  minpx = polygon[0]; maxpx = minpx;
    77   minpy = polygon[1]; maxpy = minpy; 
    78  
     77  minpy = polygon[1]; maxpy = minpy;
     78
    7979  for (i=0; i<N; i++) {
    8080    px_i = polygon[2*i];
    8181    py_i = polygon[2*i + 1];
    82  
     82
    8383    if (px_i < minpx) minpx = px_i;
    84     if (px_i > maxpx) maxpx = px_i;   
     84    if (px_i > maxpx) maxpx = px_i;
    8585    if (py_i < minpy) minpy = py_i;
    86     if (py_i > maxpy) maxpy = py_i;       
     86    if (py_i > maxpy) maxpy = py_i;
    8787  }
    88  
     88
    8989  //Begin main loop (for each point)
    9090  count = 0;
    … …  
    9393    //if (verbose){
    9494    //  if (k %((M+10)/10)==0: print 'Doing %d of %d' %(k, M)
    95        
     95
    9696    x = points[2*k];
    97     y = points[2*k + 1];       
     97    y = points[2*k + 1];
    9898
    9999    inside = 0;
    … …  
    101101    //Optimisation
    102102    if ((x > maxpx) || (x < minpx)) continue;
    103     if ((y > maxpy) || (y < minpy)) continue;       
    104 
    105     //Check polygon 
     103    if ((y > maxpy) || (y < minpy)) continue;
     104
     105    //Check polygon
    106106    for (i=0; i<N; i++) {
    107107      //printf("k,i=%d,%d\n", k, i);
    108108      j = (i+1)%N;
    109      
     109
    110110      px_i = polygon[2*i];
    111111      py_i = polygon[2*i+1];
    112112      px_j = polygon[2*j];
    113       py_j = polygon[2*j+1];     
    114      
    115       //Check for case where point is contained in line segment 
     113      py_j = polygon[2*j+1];
     114
     115      //Check for case where point is contained in line segment
    116116      if (_point_on_line(x, y, px_i, py_i, px_j, py_j)) {
    117117        if (closed == 1) {
    … …  
    120120          inside = 0;
    121121        }
    122         break; 
    123       } else {   
    124         //Check if truly inside polygon             
     122        break;
     123      } else {
     124        //Check if truly inside polygon
    125125        if ( ((py_i < y) && (py_j >= y)) ||
    126126             ((py_j < y) && (py_i >= y)) ) {
    … …  
    133133      indices[count] = k;
    134134      count++;
    135     }       
     135    }
    136136  } // End k
    137  
     137
    138138  return count;
    139139}
    … …  
    146146  // point_on_line(x, y, x0, y0, x1, y1)
    147147  //
    148  
     148
    149149  double x, y, x0, y0, x1, y1;
    150150  int res;
    151151  PyObject *result;
    152  
    153   // Convert Python arguments to C 
    154   if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1)) 
     152
     153  // Convert Python arguments to C
     154  if (!PyArg_ParseTuple(args, "dddddd", &x, &y, &x0, &y0, &x1, &y1))
    155155    return NULL;
    156156
    157  
     157
    158158  // Call underlying routine
    159159  res = _point_on_line(x, y, x0, y0, x1, y1);
    … …  
    167167
    168168PyObject *inside_polygon(PyObject *self, PyObject *args) {
    169   //def inside_polygon(point, polygon, closed, verbose, one_point): 
     169  //def inside_polygon(point, polygon, closed, verbose, one_point):
    170170  //  """Determine whether points are inside or outside a polygon
    171   // 
     171  //
    172172  //  Input:
    173173  //     point - Tuple of (x, y) coordinates, or list of tuples
    174174  //     polygon - list of vertices of polygon
    175   //     one_poin - If True Boolean value should be returned 
     175  //     one_poin - If True Boolean value should be returned
    176176  //                If False, indices of points inside returned
    177   //     closed - (optional) determine whether points on boundary should be 
    178   //     regarded as belonging to the polygon (closed = True) 
    179   //     or not (closed = False) 
    180  
    181   // 
     177  //     closed - (optional) determine whether points on boundary should be
     178  //     regarded as belonging to the polygon (closed = True)
     179  //     or not (closed = False)
     180
     181  //
    182182  //  Output:
    183183  //     If one point is considered, True or False is returned.
    184   //     If multiple points are passed in, the function returns indices 
    185   //     of those points that fall inside the polygon     
    186   //     
     184  //     If multiple points are passed in, the function returns indices
     185  //     of those points that fall inside the polygon
     186  //
    187187  //  Examples:
    188188  //     inside_polygon( [0.5, 0.5], [[0,0], [1,0], [1,1], [0,1]] )
    189   //     will evaluate to True as the point 0.5, 0.5 is inside the unit square 
    190   //     
     189  //     will evaluate to True as the point 0.5, 0.5 is inside the unit square
     190  //
    191191  //     inside_polygon( [[0.5, 0.5], [1, -0.5], [0.3, 0.2]] )
    192   //     will return the indices [0, 2] as only the first and the last point 
     192  //     will return the indices [0, 2] as only the first and the last point
    193193  //     is inside the unit square
    194   //     
     194  //
    195195  //  Remarks:
    196196  //     The vertices may be listed clockwise or counterclockwise and
    197   //     the first point may optionally be repeated.   
     197  //     the first point may optionally be repeated.
    198198  //     Polygons do not need to be convex.
    199   //     Polygons can have holes in them and points inside a hole is 
    200   //     regarded as being outside the polygon.               
    201   //           
    202   //     
    203   //  Algorithm is based on work by Darel Finley, 
    204   //  http://www.alienryderflex.com/polygon/ 
    205   //
    206   //
    207  
    208   PyArrayObject 
    209     *point, 
    210     *polygon, 
     199  //     Polygons can have holes in them and points inside a hole is
     200  //     regarded as being outside the polygon.
     201  //
     202  //
     203  //  Algorithm is based on work by Darel Finley,
     204  //  http://www.alienryderflex.com/polygon/
     205  //
     206  //
     207
     208  PyArrayObject
     209    *point,
     210    *polygon,
    211211    *indices;
    212    
    213   int closed, verbose; //Flags     
     212
     213  int closed, verbose; //Flags
    214214  int count, M, N;
    215  
    216   // Convert Python arguments to C 
    217   if (!PyArg_ParseTuple(args, "OOOii", 
    218                         &point, 
    219                         &polygon, 
    220                         &indices, 
    221                         &closed, 
    222                         &verbose)) 
     215
     216  // Convert Python arguments to C
     217  if (!PyArg_ParseTuple(args, "OOOii",
     218                        &point,
     219                        &polygon,
     220                        &indices,
     221                        &closed,
     222                        &verbose))
    223223    return NULL;
    224224
    225   M = point -> dimensions[0];    //Number of points         
    226   N = polygon -> dimensions[0];  //Number of vertices in polygon       
     225  M = point -> dimensions[0];    //Number of points
     226  N = polygon -> dimensions[0];  //Number of vertices in polygon
    227227
    228228  //printf("M=%d, N=%d\n", M, N);
    229229  // Call underlying routine
    230   count = _inside_polygon(M, N, 
     230  count = _inside_polygon(M, N,
    231231                          (double*) point -> data,
    232232                          (double*) polygon -> data,
    233233                          (int*) indices -> data,
    234                           closed, verbose); 
    235                            
     234                          closed, verbose);
     235
    236236  //NOTE: return number of points inside..
    237237  //printf("COunt=%d\n", count);
    … …  
    245245  // a,b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    246246  //
    247  
     247
    248248  double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b;
    249249  PyObject *result;
    250  
    251   // Convert Python arguments to C 
     250
     251  // Convert Python arguments to C
    252252  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
    253                         &q0, &q1, &q2)) 
     253                        &q0, &q1, &q2))
    254254    return NULL;
    255255
    256  
     256
    257257  // Call underlying routine
    258   _gradient(x0, y0, x1, y1, x2, y2, 
     258  _gradient(x0, y0, x1, y1, x2, y2,
    259259            q0, q1, q2, &a, &b);
    260260
    … …  
    271271   * three.
    272272   */
    273  
    274   //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"}, 
    275   {"gradient", gradient, METH_VARARGS, "Print out"},   
    276   {"point_on_line", point_on_line, METH_VARARGS, "Print out"},     
    277   {"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},     
     273
     274  //{"rotate", (PyCFunction)rotate, METH_VARARGS | METH_KEYWORDS, "Print out"},
     275  {"gradient", gradient, METH_VARARGS, "Print out"},
     276  {"point_on_line", point_on_line, METH_VARARGS, "Print out"},
     277  {"inside_polygon", inside_polygon, METH_VARARGS, "Print out"},
    278278  {NULL, NULL, 0, NULL}   /* sentinel */
    279279};
    … …  
    281281
    282282
    283 // Module initialisation   
     283// Module initialisation
    284284void initutil_ext(void){
    285285  Py_InitModule("util_ext", MethodTable);
    286  
    287   import_array();     //Necessary for handling of NumPY structures 
    288 }
    289 
     286
     287  import_array();     //Necessary for handling of NumPY structures
     288}
Note: See TracChangeset for help on using the changeset viewer.