source: anuga_core/source/anuga/mesh_engine/mesh_engine_c_layer.c @ 4894

Last change on this file since 4894 was 4894, checked in by duncan, 16 years ago

Good fix for the memory leak found in ticket#189

File size: 13.2 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     to return numeric arrays, check how it is done in
39     quantity_ext.c compute_gradients
40         
41     Precondition
42     End list in the pointattributelist has to have the same length
43     
44                                                 */
45
46
47#ifdef SINGLE 
48#define REAL float
49#else 
50#define REAL double
51#endif
52
53#include "triangle.h"
54
55#ifndef _STDLIB_H_
56#ifndef _WIN32
57extern "C" void *malloc(); 
58extern "C" void free(); 
59#endif
60#endif 
61
62#include "Python.h"
63
64//#define PY_ARRAY_UNIQUE_SYMBOL API_YEAH
65
66#define PY_ARRAY_UNIQUE_SYMBOL API_YEAH
67//#define NO_IMPORT_ARRAY
68#include "Numeric/arrayobject.h"
69#include <sys/types.h>
70
71 static PyObject *triang_genMesh(PyObject *self, PyObject *args){
72     
73  /* Mesh generation routine
74     pre-conditions:
75     mode must have a 'z' switch, to avoid segmentation faults
76     
77     shallow_water_ext.c gravity function
78  */
79   
80  struct triangulateio in,out;
81
82  struct triangulateio in_test;
83 
84  PyArrayObject *pointlist;
85  PyArrayObject *seglist;
86  PyArrayObject *holelist;
87  PyArrayObject *regionlist;
88  PyObject *mode;
89  PyArrayObject *pointattributelist;
90  PyArrayObject *segmarkerlist;
91  PyArrayObject *test;
92 
93 
94  PyArrayObject *r_test;
95  PyObject *R;
96
97  int dimensions[1];
98   
99  REAL Attr;
100  int i, j, iatt, n, write_here,N;
101  int a,b,c;
102  int marker;
103  int tuplesize;
104   
105  int index = 0;
106  double x,y;
107  PyObject *elems;
108  PyObject *pair;
109     
110  char *mod;
111  PyObject *holder, *holderlist, *attributelist;
112  int listsize,attsize ;
113  PyObject  *ii;
114 
115  /* used for testing numeric arrays*/
116  int n_test;     
117  if(!PyArg_ParseTuple(args,(char *)"OOOOOOOO",&pointlist,&seglist,&holelist,&regionlist,&pointattributelist,&segmarkerlist,&mode,&test)){
118    return NULL;
119  }
120 
121  /* Initialize  points */
122  in.numberofpoints =  pointlist-> dimensions[0]; 
123  in.pointlist = (double *) pointlist -> data;
124  in.pointmarkerlist = (int *)NULL; 
125 
126  /* Initialize and fill up region list */
127  in.numberofregions = regionlist-> dimensions[0]; 
128  in.regionlist = (double *) regionlist -> data;
129 
130  /*Initialize segments and segment markers */
131  in.numberofsegments = seglist -> dimensions[0]; 
132  in.segmentlist = (int *) seglist -> data;
133  /* in.segmentlist = (int *) test -> data; */
134  in.segmentmarkerlist = (int *) segmarkerlist -> data;
135 
136  /*Initialize triangles */
137  in.numberoftriangles = 0;
138  in.trianglelist = (int *)NULL;
139  in.numberoftriangleattributes = 0; 
140  in.triangleattributelist  = (REAL *)NULL; 
141     
142  /*Initialize holes */
143  in.numberofholes = holelist  -> dimensions[0];
144  if (in.numberofholes != 0) {
145    in.holelist =  (double *) holelist -> data; 
146  } else {
147    in.holelist = (REAL *)NULL;
148  }     
149 
150  /* Point attribute list */
151  if (0 == pointattributelist -> dimensions[0]) {
152    in.numberofpointattributes = 0;
153    in.pointattributelist =  NULL;
154  } else {
155    if (0 == pointattributelist -> dimensions[1]) {
156    in.numberofpointattributes = 0;
157    in.pointattributelist =  NULL;
158    } else {
159      in.numberofpointattributes = pointattributelist -> dimensions[1];
160      in.pointattributelist =  (double *) pointattributelist -> data;
161    }
162  }
163   
164  /* DEBUGGING
165  printf(" ***  hello world\n" ); 
166 
167  printf ("numberofpoints %i\n", in.numberofpoints);
168  for(i = 0; i < in.numberofpoints; i++) {
169    printf("(%f,%f)\n",in.pointlist[i* 2+ 0],in.pointlist[i* 2+ 1]);
170  }
171     
172  printf ("numberofregions %i\n", in.numberofregions);
173  for(i = 0; i < in.numberofregions; i++) {
174    printf("(%f,%f)\n",in.regionlist[i* 2+ 0],in.regionlist[i* 2+ 1]);
175  }
176   
177  printf ("numberofsegments %i\n", in.numberofsegments);
178  for(i = 0; i < in.numberofsegments; i++) {
179      printf("(%i,%i)",in.segmentlist[i* 2+ 0],in.segmentlist[i* 2+ 1]);
180      }
181 
182 
183  printf ("numberofholess %i\n", in.numberofholes);
184  for(i = 0; i < in.numberofholes; i++) {
185    printf("(%f,%f)\n",in.holelist[i* 2+ 0],in.holelist[i* 2+ 1]);
186  }
187     
188      printf ("numberoftriangles %i\n", in.numberoftriangles);
189      for(i = 0; i < in.numberoftriangles; i++) {
190      printf("(%i,%i,%i)",in.trianglelist[i* 3+ 0],in.trianglelist[i* 3+ 1], in.trianglelist[i* 3+ 2]);
191      }
192  printf(" ***  see ya world\n" );       */
193     
194  /* set up the switch arguments to the triangulation routine */
195  mod = PyString_AS_STRING(mode);
196     
197         
198  out.pointlist = (REAL *)NULL;
199  out.pointmarkerlist = (int *)NULL;
200  out.numberofpointattributes = in.numberofpointattributes;
201  out.pointattributelist = (REAL *)NULL;
202 
203  out.trianglelist = (int *)NULL;
204  out.triangleattributelist = (REAL *)NULL;                   
205  out.trianglearealist = (REAL *)NULL;                       
206  out.neighborlist = (int *)NULL;
207     
208  out.segmentlist = (int *)NULL;
209  out.segmentmarkerlist = (int *)NULL;
210     
211  out.edgelist = (int *)NULL;
212  out.edgemarkerlist = (int *)NULL;
213     
214  out.holelist = (REAL *)NULL;
215  out.regionlist = (REAL *)NULL;
216   
217 
218  /*printf("\n\nTriangulate input args: %s \n\n", mod); */
219  triangulate(mod, &in, &out, (struct triangulateio *)NULL );
220 
221 
222  /*
223  PyArray_FromDims allolws you to create a Numeric array with unitialized data.
224   The first argument is the size of the second argument (
225   the dimensions array).
226    The dimension array argument is just a 1D C array where each element of
227     the array is the size of that dimension.
228     (int dimensions[2] = { 4, 3 }; defines a 4 by 3 array.)
229     The third argument is just the desired type.
230  */
231 
232  //Py_Initialize();
233  // Testing passing a numeric array out
234  dimensions[0] = 4;
235  // Allocate space for return vectors a and b (don't DECREF)
236  r_test = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
237 
238  /* printf(" ***  back from triangulate\n" );    */
239  /*
240    ------- Pass point numbers,coordinates and neighbors back to Python ------
241    we send back a dictionary:                                               
242    { index : [ coordinates, [connections], Attribute ] }
243  */
244  holder = PyDict_New();   
245 
246  /* Add triangle list */
247  listsize = out.numberoftriangles;
248  /* printf(" out.numberoftriangles %i\n", out.numberoftriangles ); */
249  holderlist = PyList_New(listsize);
250  for(i=0; i<listsize;i++){
251    PyObject *mlist = Py_BuildValue((char *)"(i,i,i)", 
252                                    out.trianglelist[i*3  ], out.trianglelist [i*3+1], out.trianglelist [i*3+2]);   
253    PyList_SetItem(holderlist,i, mlist);
254  }   
255  ii=PyString_FromString("generatedtrianglelist");
256  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
257     
258  /* Add pointlist */
259  listsize = out.numberofpoints;
260  holderlist = PyList_New(listsize);
261     
262  for(i=0; i<out.numberofpoints;i++){
263    PyObject *mlist = Py_BuildValue((char *)"(d,d)", 
264                                    out.pointlist[i*2  ], out.pointlist[i*2+1]); 
265    PyList_SetItem(holderlist,i, mlist);
266  } 
267  ii=PyString_FromString("generatedpointlist");
268  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
269 
270  /* Add point marker list */
271  listsize = out.numberofpoints;
272  holderlist = PyList_New(listsize);
273 
274  for(i=0; i<listsize;i++){
275    PyObject *mlist = Py_BuildValue((char *)"d", 
276                                    out.pointmarkerlist[i]);     
277    PyList_SetItem(holderlist,i, mlist); 
278  } 
279  ii=PyString_FromString("generatedpointmarkerlist");
280  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
281   
282  /* Add point attribute list */
283  listsize = out.numberofpoints;
284  holderlist = PyList_New(listsize);
285 
286  for(i=0; i<listsize;i++){
287    attsize = out.numberofpointattributes;
288    attributelist = PyList_New(attsize);
289    for(iatt=0; iatt<attsize;iatt++){
290      PyObject *mlist = Py_BuildValue((char *)"d", 
291                                      out.pointattributelist[i*attsize + iatt]);     
292      PyList_SetItem(attributelist,iatt, mlist);
293    }     
294    PyList_SetItem(holderlist,i, attributelist);
295  } 
296  ii=PyString_FromString("generatedpointattributelist");
297  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist); 
298 
299 
300  /* Add triangle attribute list */
301  listsize = out.numberoftriangles;
302  holderlist = PyList_New(listsize);
303     
304  for(i=0; i<listsize; i++){
305    attsize = out.numberoftriangleattributes;
306    attributelist = PyList_New(attsize);       
307    for(iatt=0; iatt<attsize;iatt++){
308      PyObject *mlist = Py_BuildValue((char *)"d",out.triangleattributelist[i*attsize + iatt]); 
309      PyList_SetItem(attributelist,iatt, mlist);
310    }     
311    PyList_SetItem(holderlist,i, attributelist);
312  } 
313  ii=PyString_FromString("generatedtriangleattributelist");
314  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
315 
316  /* Add segment list */
317  listsize = out.numberofsegments;
318  holderlist = PyList_New(listsize);
319  for(i=0; i<listsize;i++){
320    PyObject *mlist = Py_BuildValue((char *)"(i,i)", 
321                                    out.segmentlist[i*2  ],
322                                    out.segmentlist [i*2+1]);
323    PyList_SetItem(holderlist,i, mlist);
324  }   
325  ii=PyString_FromString("generatedsegmentlist");
326  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist); 
327 
328  /* Add segment marker list */
329  listsize = out.numberofsegments;
330  holderlist = PyList_New(listsize);
331  for(i=0; i<listsize;i++){
332    PyObject *mlist = Py_BuildValue((char *)"i", 
333                                    out.segmentmarkerlist[i]);
334    PyList_SetItem(holderlist,i, mlist);
335  }   
336  ii=PyString_FromString("generatedsegmentmarkerlist");
337  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist); 
338  /* Add triangle neighbor list */
339  if (out.neighborlist != NULL) {
340    listsize = out.numberoftriangles;
341    holderlist = PyList_New(listsize);
342    for(i=0; i<listsize;i++){
343      PyObject *mlist = Py_BuildValue((char *)"(i,i,i)", 
344                                      out.neighborlist[i*3  ],
345                                      out.neighborlist [i*3+1],
346                                      out.neighborlist [i*3+2]);   
347      PyList_SetItem(holderlist,i, mlist);
348    }   
349    ii=PyString_FromString("generatedtriangleneighborlist");
350    PyDict_SetItem(holder, ii, holderlist);
351    Py_DECREF(ii); Py_DECREF(holderlist);
352  }   
353 
354 
355 
356  /* Free in/out structure memory */
357 
358  /* OUT */
359
360  if(!out.pointlist){
361    free(out.pointlist);  out.pointlist=NULL;
362  }
363  if(!out.pointmarkerlist){   
364    free(out.pointmarkerlist); out.pointmarkerlist=NULL;
365  }
366  if(!out.pointattributelist){   
367    free(out.pointattributelist); out.pointattributelist=NULL;
368  }   
369  if(!out.trianglelist){   
370    free(out.trianglelist); out.trianglelist=NULL;
371  }
372  if(!out.triangleattributelist){   
373    free(out.triangleattributelist); out.triangleattributelist=NULL;
374  }
375  if(!out.trianglearealist){   
376    free(out.trianglearealist); out.trianglearealist=NULL;
377  }
378  if(!out.neighborlist){   
379    free(out.neighborlist); out.neighborlist=NULL;
380  }
381  if(!out.segmentlist){
382    free(out.segmentlist); out.segmentlist =NULL;
383  }
384  if(!out.segmentmarkerlist){
385    free(out.segmentmarkerlist); out.segmentmarkerlist  =NULL;
386  }
387  if(!out.edgelist){
388    free(out.edgelist);  out.edgelist=NULL;
389  }
390  if(!out.edgemarkerlist){
391    free(out.edgemarkerlist);  out.edgemarkerlist=NULL;
392  }
393  if(!out.holelist){
394    free(out.holelist); out.holelist=NULL;
395  }
396  if(!out.regionlist){
397    free(out.regionlist); out.regionlist=NULL;
398  }
399 
400  /* R = Py_BuildValue((char *)"O", holder); */
401  R = Py_BuildValue((char *)"OO", holder, PyArray_Return(r_test));
402  Py_DECREF(holder); /** This fixed a  memory problem ticket#189 */
403  Py_DECREF(r_test);
404  return R;
405}
406
407/* end of my function*/
408 
409static PyMethodDef triang_methods[] = {
410  {(char *)"genMesh", triang_genMesh, 1, (char *)"Passes hull points to triangle.All arguments are lists.(<hull>,<seglist>,<holelist>,<regionlist>)-><hull>"}, 
411  {NULL,NULL}
412};
413
414void initmesh_engine_c_layer(){
415  Py_InitModule((char *)"mesh_engine_c_layer",triang_methods);
416 
417  import_array(); // Necessary for handling of NumPY structures
418}   
Note: See TracBrowser for help on using the repository browser.