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

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

working on memory freeing in the mesh generator

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