source: inundation/ga/storm_surge/pmesh/triangle/triangmodule.c @ 643

Last change on this file since 643 was 588, checked in by duncan, 20 years ago

add comments
File size: 18.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    This code is based on code by A. Pletzer
9
10    This code interfaces pmesh directly with "triangle", a general       
11    purpose triangulation code. In doing so, Python data structures         
12    are passed to C for  use by "triangle" and the results are then         
13    converted back. This was accomplished using the Python/C API.           
14                                                                           
15    I. Lists of points,edges,holes, and regions must normally be
16     created for proper triangulation. However, if there are no holes,
17     the  hole list need not be specified (however, a
18     measure of control over the shape of the convex hull is lost;
19     triangle decides).
20     
21     points list: [(x1,y1),(x2,y2),...] (Tuples of doubles)                 
22     edge list: [(Ea,Eb),(Ef,Eg),...] (Tuples of integers)                   
23     hole list: [(x1,y1),...](Tuples of doubles, one inside each hole region)
24     regionlist: [ (x1,y1,index,area),...] (Tuple of doubles)               
25     pointattributelist:[ (a1,a2..),...] (Tuple of lists)
26     mode: String - the commands sent to triangle
27         
28     This is all that is needed to accomplish an intial triangulation.       
29     A dictionary of the Triangle output is returned.                                   
30     Precondition
31     End list in the pointattributelist has to have the same length
32     
33    Compilation can be done as follows                                       
34                   
35    c++ -I/usr/local/include/python1.5 -shared -o triangmodule.so           
36    triangmodule.c triangle.o                                                */
37
38
39#ifdef SINGLE 
40#define REAL float
41#else 
42#define REAL double
43#endif 
44
45#include "triangle.h"
46
47#ifndef _STDLIB_H_
48#ifndef _WIN32
49extern "C" void *malloc(); 
50extern "C" void free(); 
51#endif
52#endif 
53
54#include "Python.h"
55
56
57 static PyObject *triang_genMesh(PyObject *self, PyObject *args){
58     
59  /* Mesh generation routine
60     pre-conditions:
61     mode must have a 'z' switch, to avoid segmentation faults
62  */
63   
64  struct triangulateio in,out;
65
66  PyObject *hull;
67  PyObject *seglist;
68  PyObject *holelist;
69  PyObject *regionlist;
70  PyObject *trianglelist;
71  PyObject *mode;
72  PyObject *firstPointAttr;
73  PyObject *pointattributelist;
74  PyObject *segmarkerlist;
75  PyObject *Attrlist;
76  REAL Attr;
77  int i, j, iatt, n;
78  int a,b,c;
79  int marker;
80  int tuplesize;
81   
82  int index = 0;
83  double x,y;
84  PyObject *elems;
85  PyObject *pair;
86     
87  char *mod;
88  PyObject *holder, *holderlist, *attributelist;
89  int listsize,attsize ;
90  PyObject  *ii;
91     
92  if(!PyArg_ParseTuple(args,(char *)"OOOOOOOO",&hull,&seglist,&holelist,&regionlist,&pointattributelist,&segmarkerlist,&trianglelist,&mode)){
93    return NULL;
94  }
95  if(!PyList_Check(hull)){
96    PyErr_SetString(PyExc_TypeError, 
97                    "incorrect first argument for 'genMesh': list required.");
98    return NULL;
99  }
100  if(!PyList_Check(hull)){
101    PyErr_SetString(PyExc_TypeError,
102                    "incorrect second argument for 'genMesh': list required.");
103    return NULL;
104  }
105  if(!PyList_Check(hull)){
106    PyErr_SetString(PyExc_TypeError,
107                    "incorrect third argument for 'genMesh': list required.");
108    return NULL;
109  }
110  if(!PyList_Check(hull)){
111    PyErr_SetString(PyExc_TypeError,
112                    "incorrect fourth argument for 'genMesh': list required.");
113    return NULL;
114  }
115  if(!PyList_Check(hull)){
116    PyErr_SetString(PyExc_TypeError,
117                    "incorrect fifth argument for 'genMesh': list required.");
118    return NULL;
119  }
120  if(!PyList_Check(hull)){
121    PyErr_SetString(PyExc_TypeError,
122                    "incorrect sixth argument for 'genMesh': list required.");
123    return NULL;
124  }
125  if(!PyList_Check(hull)){
126    PyErr_SetString(PyExc_TypeError,
127                    "incorrect seventh argument for 'genMesh': list required.");
128    return NULL;
129  }
130  if(!PyList_Check(hull)){
131    PyErr_SetString(PyExc_TypeError,
132                    "incorrect eighth argument for 'genMesh': list required.");
133    return NULL;
134  }
135     
136     
137  /* Initialize  points */
138  in.numberofpoints =  PyList_Size(hull); 
139  in.pointlist = (REAL *)malloc(in.numberofpoints*2*sizeof(REAL));
140  in.pointmarkerlist = (int *)NULL; /* malloc(in.numberofpoints*sizeof(int)); */
141     
142  /* Initialize and fill up region list */
143  in.numberofregions = PyList_Size(regionlist);
144 
145  /*printf(" numberofregions %i\n",in.numberofregions );*/
146  in.regionlist = (REAL *)malloc(in.numberofregions*4*sizeof(REAL));
147  for(i=0; i<in.numberofregions;i++){
148    elems = PyList_GetItem(regionlist,i);
149    tuplesize = PySequence_Length(elems);
150    /*printf("tuplesize %i\n",tuplesize);*/
151    if (tuplesize>0){
152      in.regionlist[4*i] = PyFloat_AsDouble(PyTuple_GetItem(elems,0));
153    } else {
154      PyErr_SetString(PyExc_TypeError,
155                      "data sent to trianglemodule error: A region has no x value.");
156      return NULL;
157    }
158    if (tuplesize>1){
159      in.regionlist[4*i+1] = PyFloat_AsDouble(PyTuple_GetItem(elems,1));
160    } else {
161      PyErr_SetString(PyExc_TypeError,
162                    "data sent to trianglemodule error: A region has no y value.");
163      return NULL;
164    }
165    if (tuplesize>2){
166      /* 2 is attribute*/
167      in.regionlist[4*i+2] = PyFloat_AsDouble(PyTuple_GetItem(elems,2));
168    } else {
169      PyErr_SetString(PyExc_TypeError,
170                    "data sent to trianglemodule error: A region has no attribute value.");
171      return NULL;
172    }
173    if (tuplesize>3){
174      in.regionlist[4*i+3] = PyFloat_AsDouble(PyTuple_GetItem(elems,3));
175    } else {
176      /* 0.0 is interpreted as no local max area by triangle  */
177      in.regionlist[4*i+3] = (REAL)0.0;
178    } 
179  }
180 
181  /*Initialize segments */
182  in.numberofsegments = PyList_Size(seglist);
183  in.segmentlist = (int *)malloc(in.numberofsegments*2*sizeof(int));
184  in.segmentmarkerlist = (int *)malloc(in.numberofsegments*sizeof(int));
185
186   
187  /*Initialize triangles */
188  in.numberofcorners = 3;
189  in.numberoftriangles = PyList_Size(trianglelist);
190  if (in.numberoftriangles  != 0) {
191    in.trianglelist = (int *)malloc(in.numberofsegments*3*sizeof(int));
192  } else {
193    in.trianglelist = (int *)NULL;
194  }
195  in.numberoftriangleattributes = 0; 
196  in.triangleattributelist  = (REAL *)NULL; 
197     
198     
199     
200  /*Initialize holes */
201  in.numberofholes = PyList_Size(holelist);
202  if (in.numberofholes != 0) {
203    in.holelist = (REAL *)malloc(in.numberofholes*2*sizeof(REAL));
204  } else {
205    in.holelist = (REAL *)NULL;
206  }     
207     
208  /* Fill up triangle list */
209  index = 0;
210  for(i = 0; i < in.numberoftriangles; i++) {
211    pair =  PyList_GetItem(trianglelist,i);
212    if (3 != PySequence_Length(pair)) {
213      PyErr_SetString(PyExc_TypeError,
214                      "data sent to trianglemodule error:Found a triangle without 3 vertices.");
215      return NULL;
216    }
217    a = (int)PyInt_AsLong(PyTuple_GetItem(pair,0));
218    in.trianglelist[index] = a;
219    index++;
220    b = (int)PyInt_AsLong(PyTuple_GetItem(pair,1));
221    in.trianglelist[index] = b;
222    index++;
223    c = (int)PyInt_AsLong(PyTuple_GetItem(pair,2));
224    in.trianglelist[index] = c;
225    index++;
226  }
227     
228     
229  /*  print the triangle list DEBUGGING
230      printf ("numberoftriangles %i\n", in.numberoftriangles);
231      for(i = 0; i < in.numberoftriangles; i++) {
232      printf("(%i,%i,%i)",in.trianglelist[i* 3+ 0],in.trianglelist[i* 3+ 1], in.trianglelist[i* 3+ 2]);
233      }*/
234     
235  /* Fill up segment list */
236  index = 0;
237  if (in.numberofsegments != PyList_Size(segmarkerlist)) {
238    PyErr_SetString(PyExc_TypeError,
239                    "data sent to trianglemodule error:number of segments doesn't = number of segment attributes.");
240    return NULL;
241  }
242   
243  for(i = 0; i < in.numberofsegments; i++) {
244    pair =  PyList_GetItem(seglist,i);
245    if (2 != PySequence_Length(pair)) {
246      PyErr_SetString(PyExc_TypeError,
247                      "data sent to trianglemodule error:Found a segment without 2 vertices.");
248      return NULL;
249    }
250    a = (int)PyInt_AsLong(PyTuple_GetItem(pair,0));
251    in.segmentlist[index] = a;
252    index++;
253    b = (int)PyInt_AsLong(PyTuple_GetItem(pair,1));
254    in.segmentlist[index] = b;
255    index++;
256    /* Fill up segment marker list */
257    marker =  (int)PyInt_AsLong(PyList_GetItem(segmarkerlist,i));
258    in.segmentmarkerlist[i] = marker;
259  }
260  /* Fill up hole list */
261  index = 0;
262  for(i = 0; i < in.numberofholes; i++) {
263    pair =  PyList_GetItem(holelist,i);
264    if (2 != PySequence_Length(pair)) {
265      PyErr_SetString(PyExc_TypeError,
266                      "data sent to trianglemodule error:Found a pair without 2 vertices.");
267      return NULL;
268    }
269    x = PyFloat_AsDouble(PyTuple_GetItem(pair,0));
270    in.holelist[index] = x;
271    index++;
272    y = PyFloat_AsDouble(PyTuple_GetItem(pair,1));
273    in.holelist[index] = y;
274    index++;
275  }
276     
277  /* Fill up point list */
278  index = 0;
279  for(i = 0; i < in.numberofpoints; i++) {
280    pair =  PyList_GetItem(hull,i);
281    if (2 != PySequence_Length(pair)) {
282      PyErr_SetString(PyExc_TypeError,
283                      "data sent to trianglemodule error:Found a vertex without an x,y values.");
284      return NULL;
285    }
286    x = PyFloat_AsDouble(PyTuple_GetItem(pair,0));
287    in.pointlist[index] = x;
288    index++;
289    y = PyFloat_AsDouble(PyTuple_GetItem(pair,1));
290    in.pointlist[index] = y;
291    index++;
292  }
293
294  /* Initialize  point attributes */
295/*      firstPointAttr = PyList_GetItem(pointattributelist,0); */
296/*      if (firstPointAttr != NULL) { */
297/*        //in.numberofpointattributes = PyList_Size(firstPointAttr); */
298/*        printf("in.numberofpointattributes =%i",PyList_Size(firstPointAttr)); */
299/*        //printf("in.numberofpointattributes =%i",in.numberofpointattributes); */
300/*      } else { */
301/*        printf("(firstPointAttr is NULL\n"); */
302/*        in.pointattributelist = NULL; */
303/*        in.numberofpointattributes = 0; */
304/*      } */
305  n = PyList_Size(pointattributelist);
306  if (n < 0) {
307    in.numberofpointattributes = 0; 
308    in.pointattributelist = (REAL *)NULL;
309  } else {
310    /* Get the # of attributes based on the size of the first point*/
311    firstPointAttr = PyList_GetItem(pointattributelist,0);
312    if (firstPointAttr == NULL) {
313      in.numberofpointattributes = 0; 
314      in.pointattributelist = (REAL *)NULL;
315    } else {
316      in.numberofpointattributes = PyList_Size(firstPointAttr);
317      in.pointattributelist = (REAL *)malloc(in.numberofpoints
318                                             *in.numberofpointattributes
319                                             *sizeof(REAL)); 
320     
321     
322      /* fill the point attribute list */
323      if (in.numberofpointattributes != 0) {
324        for(j = 0; j < in.numberofpoints; j++){
325          Attrlist = PyList_GetItem(pointattributelist,j);
326          if (in.numberofpointattributes != PySequence_Length(Attrlist)) {
327            PyErr_SetString(PyExc_TypeError,
328                            "data sent to trianglemodule error:List size of attributes is different from vertex to vertex.");
329            return NULL;
330          }
331          for(i = 0; i < in.numberofpointattributes; i++){
332            Attr = PyFloat_AsDouble(PyList_GetItem(Attrlist,i));
333            in.pointattributelist[in.numberofpointattributes*j+i]= (REAL)Attr;
334            /*   printf("i =%i\n",i);  */
335            /*   printf("j =%i\n",j);  */
336            /*   printf("in.pointattributelist[4*j+i]
337                 =%f\n",in.pointattributelist[4*j+i]); */
338          }
339        }
340      } 
341    }
342  }
343   
344     
345  /* this half from triangulate*/
346  /* set up the switch arguments to the triangulation routine */
347
348  mod = PyString_AS_STRING(mode);
349     
350         
351  out.pointlist = (REAL *)NULL;
352  out.pointmarkerlist = (int *)NULL;
353  out.numberofpointattributes = in.numberofpointattributes;
354  out.pointattributelist = (REAL *)NULL;
355 
356  out.trianglelist = (int *)NULL;
357  out.triangleattributelist = (REAL *)NULL;                   
358  out.trianglearealist = (REAL *)NULL;                       
359  out.neighborlist = (int *)NULL;
360     
361  out.segmentlist = (int *)NULL;
362  out.segmentmarkerlist = (int *)NULL;
363     
364  out.edgelist = (int *)NULL;
365  out.edgemarkerlist = (int *)NULL;
366     
367  out.holelist = (REAL *)NULL;
368  out.regionlist = (REAL *)NULL;
369   
370  /*printf("\n\nTriangulate input args: %s \n\n", mod); */
371  triangulate(mod, &in, &out, (struct triangulateio *)NULL );
372
373  /*
374    ------- Pass point numbers,coordinates and neighbors back to Python ------
375    we send back a dictionary:                                               
376    { index : [ coordinates, [connections], Attribute ] }
377  */
378  holder = PyDict_New();
379     
380  /* Add pointlist */
381  listsize = out.numberofpoints;
382  holderlist = PyList_New(listsize);
383     
384  for(i=0; i<listsize;i++){
385    PyObject *mlist = Py_BuildValue((char *)"(d,d)", 
386                                    out.pointlist[i*2  ], out.pointlist[i*2+1]); 
387    PyList_SetItem(holderlist,i, mlist);
388  } 
389  ii=PyString_FromString("generatedpointlist");
390  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
391   
392  /* Add point marker list */
393  listsize = out.numberofpoints;
394  holderlist = PyList_New(listsize);
395 
396  for(i=0; i<listsize;i++){
397    PyObject *mlist = Py_BuildValue((char *)"d", 
398                                    out.pointmarkerlist[i]);     
399    PyList_SetItem(holderlist,i, mlist); 
400  } 
401  ii=PyString_FromString("generatedpointmarkerlist");
402  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
403   
404  /* Add point attribute list */
405  listsize = out.numberofpoints;
406  holderlist = PyList_New(listsize);
407 
408  for(i=0; i<listsize;i++){
409    attsize = out.numberofpointattributes;
410    attributelist = PyList_New(attsize);
411    for(iatt=0; iatt<attsize;iatt++){
412      PyObject *mlist = Py_BuildValue((char *)"d", 
413                                      out.pointattributelist[i*attsize + iatt]);     
414      PyList_SetItem(attributelist,iatt, mlist);
415    }     
416    PyList_SetItem(holderlist,i, attributelist);
417  } 
418  ii=PyString_FromString("generatedpointattributelist");
419  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist); 
420 
421  /* Add triangle list */
422  listsize = out.numberoftriangles;
423  holderlist = PyList_New(listsize);
424  for(i=0; i<listsize;i++){
425    PyObject *mlist = Py_BuildValue((char *)"(i,i,i)", 
426                                    out.trianglelist[i*3  ], out.trianglelist [i*3+1], out.trianglelist [i*3+2]);   
427    PyList_SetItem(holderlist,i, mlist);
428  }   
429  ii=PyString_FromString("generatedtrianglelist");
430  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
431     
432     
433  /* Add triangle attribute list */
434  listsize = out.numberoftriangles;
435  holderlist = PyList_New(listsize);
436     
437  for(i=0; i<listsize; i++){
438    attsize = out.numberoftriangleattributes;
439    attributelist = PyList_New(attsize);       
440    for(iatt=0; iatt<attsize;iatt++){
441      PyObject *mlist = Py_BuildValue((char *)"d",out.triangleattributelist[i*attsize + iatt]);     
442      PyList_SetItem(attributelist,iatt, mlist);
443    }     
444    PyList_SetItem(holderlist,i, attributelist);
445  } 
446  ii=PyString_FromString("generatedtriangleattributelist");
447  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
448 
449  /* Add segment list */
450  listsize = out.numberofsegments;
451  holderlist = PyList_New(listsize);
452  for(i=0; i<listsize;i++){
453    PyObject *mlist = Py_BuildValue((char *)"(i,i)", 
454                                    out.segmentlist[i*2  ],
455                                    out.segmentlist [i*2+1]);
456    PyList_SetItem(holderlist,i, mlist);
457  }   
458  ii=PyString_FromString("generatedsegmentlist");
459  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist); 
460 
461  /* Add segment marker list */
462  listsize = out.numberofsegments;
463  holderlist = PyList_New(listsize);
464  for(i=0; i<listsize;i++){
465    PyObject *mlist = Py_BuildValue((char *)"i", 
466                                    out.segmentmarkerlist[i]);
467    PyList_SetItem(holderlist,i, mlist);
468  }   
469  ii=PyString_FromString("generatedsegmentmarkerlist");
470  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist); 
471  /* Add triangle neighbor list */
472  if (out.neighborlist != NULL) {
473    listsize = out.numberoftriangles;
474    holderlist = PyList_New(listsize);
475    for(i=0; i<listsize;i++){
476      PyObject *mlist = Py_BuildValue((char *)"(i,i,i)", 
477                                      out.neighborlist[i*3  ],
478                                      out.neighborlist [i*3+1],
479                                      out.neighborlist [i*3+2]);   
480      PyList_SetItem(holderlist,i, mlist);
481    }   
482    ii=PyString_FromString("generatedtriangleneighborlist");
483    PyDict_SetItem(holder, ii, holderlist);
484    Py_DECREF(ii); Py_DECREF(holderlist);
485  }   
486     
487  /* Free in/out structure memory */
488 
489  /* OUT */
490
491  if(!out.pointlist){
492    free(out.pointlist);  out.pointlist=NULL;
493  }
494  if(!out.pointmarkerlist){   
495    free(out.pointmarkerlist); out.pointmarkerlist=NULL;
496  }
497  if(!out.pointattributelist){   
498    free(out.pointattributelist); out.pointattributelist=NULL;
499  }
500   
501 
502  if(!out.trianglelist){   
503    free(out.trianglelist); out.trianglelist=NULL;
504  }
505  if(!out.triangleattributelist){   
506    free(out.triangleattributelist); out.triangleattributelist=NULL;
507  }
508  if(!out.trianglearealist){   
509    free(out.trianglearealist); out.trianglearealist=NULL;
510  }
511  if(!out.neighborlist){   
512    free(out.neighborlist); out.neighborlist=NULL;
513  }
514  if(!out.segmentlist){
515    free(out.segmentlist); out.segmentlist =NULL;
516  }
517  if(!out.segmentmarkerlist){
518    free(out.segmentmarkerlist); out.segmentmarkerlist  =NULL;
519  }
520  if(!out.edgelist){
521    free(out.edgelist);  out.edgelist=NULL;
522  }
523  if(!out.edgemarkerlist){
524    free(out.edgemarkerlist);  out.edgemarkerlist=NULL;
525  }
526  if(!out.holelist){
527    free(out.holelist); out.holelist=NULL;
528  }
529  if(!out.regionlist){
530    free(out.regionlist); out.regionlist=NULL;
531  }
532   
533  /*  IN */
534 
535  if(!in.pointlist){
536    free(in.pointlist);  in.pointlist   =NULL;             
537  }
538  if(!in.pointattributelist){
539    free(in.pointattributelist); in.pointattributelist =NULL;       
540  }
541  if(!in.pointmarkerlist){
542    free(in.pointmarkerlist); in.pointmarkerlist    =NULL;       
543  }
544  if(!in.segmentlist){
545    free(in.segmentlist);  in.segmentlist    =NULL;         
546  }
547  if(!in.segmentmarkerlist){
548    free(in.segmentmarkerlist); in.segmentmarkerlist    =NULL;   
549  }
550  if(!in.regionlist){
551    free(in.regionlist); in.regionlist  =NULL;           
552  }
553  if(!in.holelist){
554    free(in.holelist); in.holelist=NULL;
555  }
556
557  //return Py_BuildValue((char *)"O", holder);
558  return Py_BuildValue((char *)"O", holder);
559}
560
561/* end of my function*/
562
563
564   
565   
566static PyMethodDef triang_methods[] = {
567  {(char *)"genMesh", triang_genMesh, 1, (char *)"Passes hull points to triangle.All arguments are lists.(<hull>,<seglist>,<holelist>,<regionlist>)-><hull>"}, 
568  {NULL,NULL}
569};
570
571void inittriang(){
572  Py_InitModule((char *)"triang",triang_methods);
573}   
Note: See TracBrowser for help on using the repository browser.