source: branches/numpy/anuga/mesh_engine/mesh_engine_c_layer.c @ 6410

Last change on this file since 6410 was 6410, checked in by rwilson, 15 years ago

numpy changes.

File size: 12.7 KB
Line 
1/* this appears to be necessary to run on Linux */
2
3#undef _STDLIB_H_
4#define _STDLIB_H_
5
6
7/*
8
9    This code interfaces pmesh directly with "triangle", a general       
10    purpose triangulation code. In doing so, Python numeric data structures
11     are passed to C for  use by "triangle" and the results are then         
12    converted back. This was accomplished using the Python/C API.           
13                                                                           
14    I. Arrays of points,edges,holes, and regions must normally be
15     created for proper triangulation. However, if there are no holes,
16     the  hole list need not be specified (however, a
17     measure of control over the shape of the convex hull is lost;
18     triangle decides).
19     
20     points array: [(x1,y1),(x2,y2),...] ( doubles)                 
21     edge array: [(Ea,Eb),(Ef,Eg),...] ( integers)                   
22     hole array: [(x1,y1),...]( doubles, one inside each hole region)
23     region array: [ (x1,y1,index,area),...] ( doubles)               
24     pointattribute array:[ (a1,a2..),...] (doubles)
25     mode: String - the commands sent to triangle
26         
27     This is all that is needed to accomplish an intial triangulation.       
28     A dictionary of the Triangle output is returned.
29           
30     I thought this function was leaking. That the returned data
31     structure didn't get garbage collected.  I decided this using 
32     gc.get_objects() and seeing the dic still hanging around.
33     I'm not so sure this means its leaking.  Check by looking at
34     the overall memory use.  Anyhow, I delete the lists that are
35     returned in the dic structre after they are used now, in alpha
36     shape and mesh_engine.
37
38         
39     Precondition
40     End list in the pointattributelist has to have the same length
41     
42                                                 */
43
44
45#ifdef SINGLE 
46#define REAL float
47#else 
48#define REAL double
49#endif
50
51#include "triangle.h"
52
53#ifndef _STDLIB_H_
54#ifndef _WIN32
55extern "C" void *malloc(); 
56extern "C" void free(); 
57#endif
58#endif 
59
60#include "Python.h"
61
62#include "util_ext.h"
63
64
65//#define PY_ARRAY_UNIQUE_SYMBOL API_YEAH
66
67#define PY_ARRAY_UNIQUE_SYMBOL API_YEAH
68//#define NO_IMPORT_ARRAY
69#include "numpy/arrayobject.h"
70#include <sys/types.h>
71
72 static PyObject *triang_genMesh(PyObject *self, PyObject *args){
73     
74  /* Mesh generation routine
75     pre-conditions:
76     mode must have a 'z' switch, to avoid segmentation faults
77     
78     shallow_water_ext.c gravity function
79  */
80   
81  struct triangulateio in,out;
82
83  struct triangulateio in_test;
84 
85  PyArrayObject *pointlist;
86  PyArrayObject *seglist;
87  PyArrayObject *holelist;
88  PyArrayObject *regionlist;
89  PyObject *mode;
90  PyArrayObject *pointattributelist;
91  PyArrayObject *segmarkerlist;
92  PyArrayObject *test;
93 
94 
95  PyArrayObject *gentrianglelist;
96  PyArrayObject *genpointlist;
97  PyArrayObject *genseglist;
98  PyArrayObject *genpointmarkerlist;
99  PyArrayObject *genpointattributelist;
100  PyArrayObject *gentriangleattributelist;
101  PyArrayObject *gensegmentlist;
102  PyArrayObject *gensegmentmarkerlist;
103  PyArrayObject *genneighborlist;
104  PyObject *R;
105
106  int dimensions[2];
107   
108  REAL Attr;
109  int i, j, iatt, n, write_here,N;
110  int a,b,c;
111  int marker;
112  int tuplesize;
113   
114  int index = 0;
115  double x,y;
116  PyObject *elems;
117  PyObject *pair;
118     
119  char *mod;
120   
121  if(!PyArg_ParseTuple(args,(char *)"OOOOOOO",&pointlist,&seglist,&holelist,&regionlist,&pointattributelist,&segmarkerlist,&mode)){
122    return NULL;
123  }
124
125  // check that numpy array objects arrays are C contiguous memory
126  CHECK_C_CONTIG(pointlist);
127  CHECK_C_CONTIG(seglist);
128  CHECK_C_CONTIG(holelist);
129  CHECK_C_CONTIG(regionlist);
130  CHECK_C_CONTIG(pointattributelist);
131  CHECK_C_CONTIG(segmarkerlist);
132 
133  /* Initialize  points */
134  in.numberofpoints =  pointlist-> dimensions[0]; 
135  in.pointlist = (double *) pointlist -> data;
136  in.pointmarkerlist = (int *)NULL; 
137 
138  /* Initialize and fill up region list */
139  in.numberofregions = regionlist-> dimensions[0]; 
140 
141  in.regionlist = (double *) regionlist -> data;
142 
143  /*Initialize segments and segment markers */
144  in.numberofsegments = seglist -> dimensions[0]; 
145  in.segmentlist = (int *) seglist -> data;
146  in.segmentmarkerlist = (int *) segmarkerlist -> data;
147  //in.segmentmarkerlist = (int *) NULL;
148 
149  /*Initialize triangles */
150  in.numberoftriangles = 0;
151  in.trianglelist = (int *)NULL;
152  in.numberoftriangleattributes = 0; 
153  in.triangleattributelist  = (REAL *)NULL; 
154  in.trianglearealist  = (REAL *)NULL; 
155  in.neighborlist = (int *)NULL;
156  in.numberofcorners = 0;
157     
158  /*Initialize holes */
159  in.numberofholes = holelist  -> dimensions[0];
160  if (in.numberofholes != 0) {
161    in.holelist =  (double *) holelist -> data; 
162  } else {
163    in.holelist = (REAL *)NULL;
164  }     
165 
166  /* Point attribute list */
167   /*printf ("in.pointattributelist -> dimensions[0] %i\n", pointattributelist -> dimensions[0]);
168  printf ("in.pointattributelist -> dimensions[1] %i\n", pointattributelist -> dimensions[1]); */
169  if (0 == pointattributelist -> dimensions[0]) {
170    in.numberofpointattributes = 0;
171    in.pointattributelist =  (double *) NULL;
172  } else {
173    if (0 == pointattributelist -> dimensions[1]) {
174    in.numberofpointattributes = 0;
175    in.pointattributelist =  (double *) NULL;
176    } else {
177      in.numberofpointattributes = pointattributelist -> dimensions[1];
178      in.pointattributelist =  (double *) pointattributelist -> data;
179    }
180  }
181   
182  /* DEBUGGING
183  printf(" ***  hello world\n" ); 
184 
185  printf ("numberofpoints %i\n", in.numberofpoints);
186  printf ("numberofpointattributes %i\n", in.numberofpointattributes);
187  for(i = 0; i < in.numberofpoints; i++) {
188    printf("(%f,%f)\n",in.pointlist[i* 2+ 0],in.pointlist[i* 2+ 1]);
189    for(j = 0; j < in.numberofpointattributes; j++) {
190      printf("point att (%f)\n",in.pointattributelist[i* in.numberofpointattributes + j]);
191    }
192   
193  }
194     
195  printf ("numberofregions %i\n", in.numberofregions);
196  for(i = 0; i < in.numberofregions; i++) {
197    printf("(%f,%f)  ",in.regionlist[i* 4+ 0],in.regionlist[i* 4+ 1]);
198    printf("index %f Area %f\n",in.regionlist[i* 4+ 2],in.regionlist[i* 4+ 3]);
199  }
200   
201  printf ("numberofsegments %i\n", in.numberofsegments);
202  for(i = 0; i < in.numberofsegments; i++) {
203      printf("(%i,%i)\n",in.segmentlist[i* 2+ 0],in.segmentlist[i* 2+ 1]);
204      printf("Segment marker (%i)\n",in.segmentmarkerlist[i + 0]);
205      }
206 
207 
208  printf ("numberofholess %i\n", in.numberofholes);
209  for(i = 0; i < in.numberofholes; i++) {
210    printf("(%f,%f)\n",in.holelist[i* 2+ 0],in.holelist[i* 2+ 1]);
211  }
212     
213      printf ("numberoftriangles %i\n", in.numberoftriangles);
214      for(i = 0; i < in.numberoftriangles; i++) {
215      printf("(%i,%i,%i)",in.trianglelist[i* 3+ 0],in.trianglelist[i* 3+ 1], in.trianglelist[i* 3+ 2]);
216      }
217  printf(" ***  see ya world\n" );       
218      */
219  /* set up the switch arguments to the triangulation routine */
220  mod = PyString_AS_STRING(mode);
221     
222         
223  out.pointlist = (REAL *)NULL;
224  out.pointmarkerlist = (int *)NULL;
225  //out.numberofpointattributes = in.numberofpointattributes;
226  out.pointattributelist = (REAL *)NULL;
227 
228  out.trianglelist = (int *)NULL;
229  out.triangleattributelist = (REAL *)NULL;                   
230  out.trianglearealist = (REAL *)NULL;                       
231  out.neighborlist = (int *)NULL;
232     
233  out.segmentlist = (int *)NULL;
234  out.segmentmarkerlist = (int *)NULL;
235     
236  out.edgelist = (int *)NULL;
237  out.edgemarkerlist = (int *)NULL;
238     
239  out.holelist = (REAL *)NULL;
240  out.regionlist = (REAL *)NULL;
241   
242 
243  triangulate(mod, &in, &out, (struct triangulateio *)NULL );
244 
245 
246   
247  /*
248    ------- Pass point numbers,coordinates and neighbors back to Python ------
249    we send back a dictionary:                                               
250    { index : [ coordinates, [connections], Attribute ] }
251  */ 
252 
253  /*
254 
255     to return numeric arrays, check how it is done in
256     abs quantity_ext.c compute_gradients
257     
258  PyArray_FromDims allolws you to create a numeric array with unitialized data.
259   The first argument is the size of the second argument (
260   the dimensions array).
261    The dimension array argument is just a 1D C array where each element of
262     the array is the size of that dimension.
263     (int dimensions[2] = { 4, 3 }; defines a 4 by 3 array.)
264     The third argument is just the desired type.
265  */
266 
267 
268  /* Add triangle list */
269  dimensions[0] = out.numberoftriangles;
270  dimensions[1] = 3;   
271  gentrianglelist = (PyArrayObject *) PyArray_FromDimsAndData(2,
272                                            dimensions,
273                                            PyArray_INT, 
274                                            (char*) out.trianglelist); 
275   
276  /* Add pointlist */
277  dimensions[0] = out.numberofpoints;
278  dimensions[1] = 2;   
279  genpointlist = (PyArrayObject *) PyArray_FromDimsAndData(2,
280                                         dimensions,
281                                         PyArray_DOUBLE, 
282                                         (char*) out.pointlist);
283                                           
284 
285  /* Add point marker list */
286  dimensions[0] = out.numberofpoints;
287  genpointmarkerlist = (PyArrayObject *) PyArray_FromDimsAndData(1, 
288                                         dimensions, 
289                                         PyArray_INT,
290                                        (char*) out.pointmarkerlist);
291 
292  /* Add point attribute list */
293  dimensions[0] = out.numberofpoints;
294  dimensions[1] = out.numberofpointattributes;   
295  genpointattributelist = (PyArrayObject *) PyArray_FromDimsAndData(2, 
296                                          dimensions, 
297                                          PyArray_DOUBLE,
298                                          (char*) out.pointattributelist);
299 
300 
301 
302  /* Add triangle attribute list */
303  dimensions[0] = out.numberoftriangles;
304  dimensions[1] = out.numberoftriangleattributes;   
305  gentriangleattributelist = (PyArrayObject *) PyArray_FromDimsAndData(2, 
306                                           dimensions, 
307                                           PyArray_DOUBLE,
308                                          (char*)out.triangleattributelist);
309 
310  /* Add segment list */
311  dimensions[0] = out.numberofsegments;
312  dimensions[1] = 2;   
313  gensegmentlist = (PyArrayObject *) PyArray_FromDimsAndData(2, 
314                                                    dimensions, 
315                                                    PyArray_INT,
316                                                    (char*)out.segmentlist);
317 
318 
319  /* Add segment marker list */
320  dimensions[0] = out.numberofsegments;
321  gensegmentmarkerlist = (PyArrayObject *) PyArray_FromDimsAndData(1, 
322                                            dimensions, 
323                                            PyArray_INT,
324                                           (char*)out.segmentmarkerlist);
325 
326  /* Add triangle neighbor list */
327  if (out.neighborlist != NULL) {
328    dimensions[0] = out.numberoftriangles;
329    dimensions[1] = 3;   
330    genneighborlist = (PyArrayObject *) PyArray_FromDimsAndData(2, 
331                                                 dimensions, 
332                                                 PyArray_INT,
333                                                 (char*)out.neighborlist);
334  }else{ 
335    dimensions[0] = 0;
336    dimensions[1] = 0;   
337    genneighborlist = (PyArrayObject *) PyArray_FromDims(2, 
338                                                         dimensions, 
339                                                         PyArray_INT);
340  }
341 
342 
343  R = Py_BuildValue((char *)"OOOOOOOO"
344                    ,PyArray_Return(gentrianglelist)
345                    ,PyArray_Return(genpointlist)
346                    ,PyArray_Return(genpointmarkerlist)
347                    ,PyArray_Return(genpointattributelist)
348                    ,PyArray_Return(gentriangleattributelist)
349                    ,PyArray_Return(gensegmentlist)
350                    ,PyArray_Return(gensegmentmarkerlist)
351                    ,PyArray_Return(genneighborlist)
352                    );
353                   
354  Py_DECREF(gentrianglelist);
355  Py_DECREF(genpointlist);
356  Py_DECREF(genpointmarkerlist);
357  Py_DECREF(genpointattributelist);
358  Py_DECREF(gentriangleattributelist);
359  Py_DECREF(gensegmentlist);
360  Py_DECREF(gensegmentmarkerlist);
361  Py_DECREF(genneighborlist);
362 
363 
364  /* These memory blocks are passed into numeric arrays
365   so don't free them  */
366  /*
367  if(!out.trianglelist){   
368    free(out.trianglelist); out.trianglelist=NULL;
369   
370   
371  if(!out.pointlist){
372    free(out.pointlist);  out.pointlist=NULL;
373  } 
374 
375  if(!out.pointmarkerlist){   
376    free(out.pointmarkerlist); out.pointmarkerlist=NULL;
377  }
378 
379  if(!out.pointattributelist){   
380    free(out.pointattributelist); out.pointattributelist=NULL;
381  }   
382 
383  if(!out.triangleattributelist){   
384    free(out.triangleattributelist); out.triangleattributelist=NULL;
385  }
386  if(!out.segmentlist){
387    free(out.segmentlist); out.segmentlist =NULL;
388  }
389  if(!out.segmentmarkerlist){
390    free(out.segmentmarkerlist); out.segmentmarkerlist  =NULL;
391  }
392  if(!out.neighborlist){   
393    free(out.neighborlist); out.neighborlist=NULL;
394  }
395   */
396
397   
398 
399  /* Free in/out structure memory */
400 
401  if(!out.trianglearealist){   
402    free(out.trianglearealist); out.trianglearealist=NULL;
403  }
404  if(!out.edgelist){
405    free(out.edgelist);  out.edgelist=NULL;
406  }
407  if(!out.edgemarkerlist){
408    free(out.edgemarkerlist);  out.edgemarkerlist=NULL;
409  }
410  if(!out.holelist){
411    free(out.holelist); out.holelist=NULL;
412  }
413  if(!out.regionlist){
414    free(out.regionlist); out.regionlist=NULL;
415  }
416  return R;
417}
418
419/* end of my function*/
420 
421static PyMethodDef triang_methods[] = {
422  {(char *)"genMesh", triang_genMesh, 1, (char *)"Passes hull points to triangle.All arguments are lists.(<hull>,<seglist>,<holelist>,<regionlist>)-><hull>"}, 
423  {NULL,NULL}
424};
425
426void initmesh_engine_c_layer(){
427  Py_InitModule((char *)"mesh_engine_c_layer",triang_methods);
428 
429  import_array(); // Necessary for handling of NumPY structures
430}   
Note: See TracBrowser for help on using the repository browser.