Changeset 8504


Ignore:
Timestamp:
Aug 11, 2012, 10:48:40 PM (13 years ago)
Author:
steve
Message:

Moved some slow code to C in buid_neighbour_structure

Location:
trunk/anuga_core/source/anuga/abstract_2d_finite_volumes
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/hashtable.c

    r8503 r8504  
    2020//==============================================================================
    2121
    22 struct edge {
    23     int key;                    /* key of form i m + j */
     22typedef struct  {
     23    int i;
     24    int j;
     25} edge_key_t;
     26
     27typedef struct {
     28    edge_key_t key;              /* key of form i , j */
    2429    int vol_id;                /* id of vol containing this edge */
    2530    int edge_id;               /* edge_id of edge in this vol */
    2631    UT_hash_handle hh;         /* makes this structure hashable */
    27 };
    28 
    29 struct edge *edgetable = NULL;
    30 
    31 void add_edge(int edge_key, int vol_id, int edge_id) {
    32     struct edge *s;
    33 
    34     s = (struct edge*)malloc(sizeof(struct edge));
    35     s->key = edge_key;
     32} edge_t;
     33
     34edge_t *edgetable = NULL;
     35
     36void add_edge(edge_key_t key, int vol_id, int edge_id) {
     37    edge_t *s;
     38
     39    s = (edge_t*) malloc(sizeof(edge_t));
     40    memset(s, 0, sizeof(edge_t));
     41    s->key.i = key.i;
     42    s->key.j = key.j;
    3643    s->vol_id = vol_id;
    3744    s->edge_id = edge_id;
    38     HASH_ADD_INT( edgetable, key, s );  /* key: name of key field */
    39 }
    40 
    41 struct edge *find_edge(int edge_key) {
    42     struct edge *s;
    43 
    44     HASH_FIND_INT( edgetable, &edge_key, s );  /* s: output pointer */
     45    HASH_ADD(hh, edgetable, key, sizeof(edge_key_t), s);  /* key: name of key field */
     46}
     47
     48edge_t *find_edge(edge_key_t key) {
     49    edge_t *s;
     50
     51    HASH_FIND(hh, edgetable, &key, sizeof(edge_key_t), s);  /* s: output pointer */
    4552    return s;
    4653}
    4754
    48 void delete_edge(struct edge *edge) {
     55void delete_edge(edge_t *edge) {
    4956    HASH_DEL( edgetable, edge);  /* user: pointer to deletee */
    5057    free(edge);
    … …  
    5259
    5360void delete_all() {
    54   struct edge *current_edge, *tmp;
     61  edge_t *current_edge, *tmp;
    5562
    5663  HASH_ITER(hh, edgetable, current_edge, tmp) {
    57     HASH_DEL(edgetable,current_edge);  /* delete it (edgetable advances to next) */
     64    HASH_DEL(edgetable, current_edge);  /* delete it (edgetable advances to next) */
    5865    free(current_edge);            /* free it */
    5966  }
    … …  
    6168
    6269void print_edges() {
    63     struct edge *s;
    64 
    65     for(s=edgetable; s != NULL; s=(struct edge*)(s->hh.next)) {
    66         printf("edge key %d: vol_id %d  edge_id %d\n", s->key, s->vol_id, s->edge_id);
     70    edge_t *s;
     71
     72    for(s=edgetable; s != NULL; s=(edge_t*)(s->hh.next)) {
     73        printf("edge key i %d i %d vol_id %d  edge_id %d\n",
     74                      s->key.i, s->key.j, s->vol_id, s->edge_id);
    6775    }
    6876}
    6977
    70 int vol_id_sort(struct edge *a, struct edge *b) {
     78int vol_id_sort(edge_t *a, edge_t *b) {
    7179    return (a->vol_id - b->vol_id);
    7280}
    7381
    74 int key_sort(struct edge *a, struct edge *b) {
    75     return (a->key - b->key);
     82int key_sort(edge_t *a, edge_t *b) {
     83    return (a->key.i - b->key.i);
    7684}
    7785
    … …  
    8896// Code to calculate neighbour structure
    8997//==============================================================================
    90 int _create_neighbours(int N,
     98int _build_neighbour_structure(int N, int M,
     99                      long* triangles,
    91100                      long* neighbours,
    92101                      long* neighbour_edges,
    93                       long* number_of_boundaries,
    94                       long* triangles) {
     102                      long* number_of_boundaries)
     103                      {
    95104    int k;
    96105    int k3;
    97106    int n0,n1,n2;
    98     int key;
    99107    int vol_id;
    100108    int edge_id;
    101     struct edge *s;
     109    int err = 0;
     110    edge_t *s;
     111    edge_key_t key;
    102112
    103113    //--------------------------------------------------------------------------
    … …  
    106116    // two nodes defining the edge
    107117    //--------------------------------------------------------------------------
    108     for (k=0; k<N; k++) {
     118    for (k=0; k<M; k++) {
    109119
    110120        // Run through triangles
    … …  
    117127        // Add edges
    118128
     129        //----------------------------------------------------------------------
    119130        // edge 0, n1 - n2
    120         key = n1*N + n2;
     131        //----------------------------------------------------------------------
     132        key.i = n1;
     133        key.j = n2;
    121134        vol_id = k;
    122135        edge_id = 0;
    123 
     136       
     137        // Error if duplicates
     138        s = find_edge(key);
     139        if (s) {
     140            err = 1;
     141            break;
     142        }
     143        // Otherwise add edge
    124144        add_edge(key, vol_id, edge_id);
    125145
     146        //----------------------------------------------------------------------
    126147        // edge 1, n2 - n0
    127         key = n2*N + n0;
     148        //----------------------------------------------------------------------
     149        key.i = n2;
     150        key.j = n0;
    128151        vol_id = k;
    129152        edge_id = 1;
    130153
     154        // Error if duplicates
     155        s = find_edge(key);
     156        if (s) {
     157            err = 1;
     158            break;
     159        }
     160        // Otherwise add edge
    131161        add_edge(key, vol_id, edge_id);
    132162
     163        //----------------------------------------------------------------------
    133164        // edge 2, n0 - n1
    134         key = n0*N + n1;
     165        //----------------------------------------------------------------------
     166        key.i = n0;
     167        key.j = n1;
    135168        vol_id = k;
    136169        edge_id = 2;
    137170
     171        // Error if duplicates
     172        s = find_edge(key);
     173        if (s) {
     174            err = 1;
     175            break;
     176        }
     177        // Otherwise add edge
    138178        add_edge(key, vol_id, edge_id);
    139179
    140180    }
     181
     182   
     183    //--------------------------------------------------------------------------
     184    // return with an error code if duplicate key found
     185    // Clean up hashtable
     186    //--------------------------------------------------------------------------
     187    if (err) {
     188        //printf("Duplicate Keys:\n");
     189        //printf("key.i %d key.j %d vol_id %d edge_id %d \n",
     190        //       s->key.i, s->key.j,s->vol_id,s->edge_id);
     191        //printf("key.i %d key.j %d vol_id %d edge_id %d \n",
     192        //       key.i,key.j,vol_id,edge_id);
     193        delete_all();
     194        return err;
     195    }
     196
    141197
    142198    //--------------------------------------------------------------------------
    … …  
    145201    //reverse direction of segments and lookup neighbours.
    146202    //--------------------------------------------------------------------------
    147     for (k=0; k<N; k++) {
     203    for (k=0; k<M; k++) {
    148204
    149205        // Run through triangles
    … …  
    159215
    160216        // edge 0, n1 - n2
    161         key = n2*N + n1;
     217        key.i = n2;
     218        key.j = n1;
     219        s = find_edge(key);
     220        if (s) {
     221            neighbours[k3]      = s -> vol_id;
     222            neighbour_edges[k3] = s -> edge_id;
     223            number_of_boundaries[k] -= 1;
     224        }
     225
     226        // edge 1, n2 - n0
     227        key.i = n0;
     228        key.j = n2;
     229        s = find_edge(key);
     230        if (s) {
     231            neighbours[k3+1]      = s -> vol_id;
     232            neighbour_edges[k3+1] = s -> edge_id;
     233            number_of_boundaries[k] -= 1;
     234        }
     235
     236        // edge 2, n0 - n1
     237        key.i = n1;
     238        key.j = n0;
    162239        s = find_edge(key);
    163240        if (s) {
    164241            neighbours[k3+2]      = s -> vol_id;
    165242            neighbour_edges[k3+2] = s -> edge_id;
    166             number_of_boundaries[k3+2] -= 1;
    167         }
    168 
    169         // edge 1, n2 - n0
    170         key = n0*N + n2;
    171         s = find_edge(key);
    172         if (s) {
    173             neighbours[k3+2]      = s -> vol_id;
    174             neighbour_edges[k3+2] = s -> edge_id;
    175             number_of_boundaries[k3+2] -= 1;
    176         }
    177 
    178         // edge 2, n0 - n1
    179         key = n1*N + n0;
    180         s = find_edge(key);
    181         if (s) {
    182             neighbours[k3+2]      = s -> vol_id;
    183             neighbour_edges[k3+2] = s -> edge_id;
    184             number_of_boundaries[k3+2] -= 1;
     243            number_of_boundaries[k] -= 1;
    185244        }
    186245
    … …  
    189248    delete_all();  /* free any structures */
    190249
    191     return 0;
     250    return err;
    192251
    193252}
    … …  
    198257//==============================================================================
    199258
    200 PyObject *create_neighbours(PyObject *self, PyObject *args) {
     259PyObject *build_neighbour_structure(PyObject *self, PyObject *args) {
    201260
    202261  /*
    203262   * Update neighbours array using triangles array
     263   *
     264   * N is number of nodes (vertices)
     265   * triangle nodes defining triangles
     266   * neighbour across edge_id
     267   * neighbour_edges edge_id of edge in neighbouring triangle
     268   * number_of_boundaries
    204269  */
    205270
    … …  
    207272        PyArrayObject *triangles;
    208273
    209         int N, err;
     274        int N; // Number of nodes (read in)
     275        int M; // Number of triangles (calculated from triangle array)
     276        int err;
    210277
    211278
    212279        // Convert Python arguments to C
    213         if (!PyArg_ParseTuple(args, "OOOO", &neighbours,
     280        if (!PyArg_ParseTuple(args, "iOOOO", &N,
     281                                            &triangles,
     282                                            &neighbours,
    214283                                            &neighbour_edges,
    215                                             &number_of_boundaries,
    216                                             &triangles)) {
     284                                            &number_of_boundaries
     285                                            )) {
    217286          PyErr_SetString(PyExc_RuntimeError,
    218287                          "hashtable.c: create_neighbours could not parse input");
    … …  
    220289        }
    221290
     291        CHECK_C_CONTIG(triangles);
    222292        CHECK_C_CONTIG(neighbours);
    223293        CHECK_C_CONTIG(neighbour_edges);
    224294        CHECK_C_CONTIG(number_of_boundaries);
    225         CHECK_C_CONTIG(triangles);
    226 
    227         N = triangles -> dimensions[0];
    228 
    229         err = _create_neighbours(N,
     295
     296
     297        M = triangles -> dimensions[0];
     298
     299
     300        err = _build_neighbour_structure(N, M,
     301                      (long*) triangles  -> data,
    230302                      (long*) neighbours -> data,
    231303                      (long*) neighbour_edges -> data,
    232                       (long*) number_of_boundaries -> data,
    233                       (long*) triangles  -> data);
     304                      (long*) number_of_boundaries -> data);
    234305
    235306
    236307        if (err != 0) {
    237308          PyErr_SetString(PyExc_RuntimeError,
    238                           "hashtable.c: error creating neighbours");
     309                          "Duplicate Edge");
    239310          return NULL;
    240311        }
    241312
    242313        // Release and return
    243         Py_DECREF(triangles);
    244         Py_DECREF(neighbours);
    245         Py_DECREF(neighbour_edges);
    246         Py_DECREF(number_of_boundaries);
     314        //Py_DECREF(triangles);
     315        //Py_DECREF(neighbours);
     316        //Py_DECREF(neighbour_edges);
     317        //Py_DECREF(number_of_boundaries);
    247318
    248319        return Py_BuildValue("");
    … …  
    256327// Method table for python module
    257328static struct PyMethodDef MethodTable[] = {
    258         {"create_neighbours", create_neighbours, METH_VARARGS, "Print out"},
     329        {"build_neighbour_structure", build_neighbour_structure, METH_VARARGS, "Print out"},
    259330        {NULL, NULL, 0, NULL}   // sentinel
    260331};

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r8503 r8504  
    233233
    234234
    235     def build_neighbour_structure(self):
     235    def build_neighbour_structure_python(self):
    236236        """Update all registered triangles to point to their neighbours.
    237237
    … …  
    256256            b = self.triangles[i, 1]
    257257            c = self.triangles[i, 2]
    258 #            if neighbourdict.has_key((a,b)):
    259 #                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
    260 #                    raise Exception(msg)
    261 #            if neighbourdict.has_key((b,c)):
    262 #                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
    263 #                    raise Exception(msg)
    264 #            if neighbourdict.has_key((c,a)):
    265 #                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
    266 #                    raise Exception(msg)
     258            if neighbourdict.has_key((a,b)):
     259                    msg = "Edge 2 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[a,b][1],neighbourdict[a,b][0])
     260                    raise Exception(msg)
     261            if neighbourdict.has_key((b,c)):
     262                    msg = "Edge 0 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[b,c][1],neighbourdict[b,c][0])
     263                    raise Exception(msg)
     264            if neighbourdict.has_key((c,a)):
     265                    msg = "Edge 1 of triangle %d is duplicating edge %d of triangle %d.\n" %(i,neighbourdict[c,a][1],neighbourdict[c,a][0])
     266                    raise Exception(msg)
    267267
    268268            neighbourdict[a,b] = (i, 2) #(id, edge)
    … …  
    295295                self.number_of_boundaries[i] -= 1
    296296
     297    def build_neighbour_structure(self):
     298        """Update all registered triangles to point to their neighbours.
     299
     300        Also, keep a tally of the number of boundaries for each triangle
     301
     302        Postconditions:
     303          neighbours and neighbour_edges is populated
     304          number_of_boundaries integer array is defined.
     305        """
     306
     307        import hashtable
     308       
     309        N = self.number_of_nodes
     310
     311       
     312        hashtable.build_neighbour_structure(N,
     313                                            self.triangles,
     314                                            self.neighbours,
     315                                            self.neighbour_edges,
     316                                            self.number_of_boundaries)
     317 
    297318
    298319    def build_surrogate_neighbour_structure(self):
    … …  
    312333
    313334        N = len(self) #Number of triangles
    314         for i in range(N):
     335        for i in xrange(N):
    315336            #Find all neighbouring volumes that are not boundaries
    316             for k in range(3):
     337            for k in xrange(3):
    317338                if self.neighbours[i, k] < 0:
    318339                    self.surrogate_neighbours[i, k] = i #Point this triangle
    … …  
    336357        if boundary is None:
    337358            boundary = {}
    338             for vol_id in range(len(self)):
    339                 for edge_id in range(0, 3):
     359            for vol_id in xrange(len(self)):
     360                for edge_id in xrange(0, 3):
    340361                    if self.neighbours[vol_id, edge_id] < 0:
    341362                        boundary[(vol_id, edge_id)] = default_boundary_tag
    … …  
    352373
    353374            #Check that all boundary segments are assigned a tag
    354             for vol_id in range(len(self)):
    355                 for edge_id in range(0, 3):
     375            for vol_id in xrange(len(self)):
     376                for edge_id in xrange(0, 3):
    356377                    if self.neighbours[vol_id, edge_id] < 0:
    357378                        if not boundary.has_key( (vol_id, edge_id) ):
Note: See TracChangeset for help on using the changeset viewer.