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

Last change on this file since 541 was 349, checked in by duncan, 21 years ago

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