source: issm/trunk/src/mex/TriMesh/TriMesh.cpp@ 1904

Last change on this file since 1904 was 1904, checked in by Eric.Larour, 16 years ago

Valgrind memory checker now working.
File size: 10.3 KB
Line 
1/*
2 * TriMesh: out of a domain outline file ( Argus format ),
3 * use the Triangle package to create a triangular mesh
4 *
5 */
6
7#include "./TriMesh.h"
8
9
10void mexFunction( int nlhs, mxArray* plhs[],
11 int nrhs, const mxArray* prhs[] )
12{
13
14
15 /*Matlab arrays: */
16 mxArray* pmxa_array=NULL;
17 int i,j;
18 int counter,counter2,backcounter;
19 int prhs_counter;
20
21 /* returned quantities: */
22
23 double* index=NULL;
24 double* x=NULL;
25 double* y=NULL;
26 double* segments=NULL;
27 double* segmentmarkerlist=NULL;
28
29 /* input: */
30 char* domainname=NULL;
31 char* riftname=NULL;
32 double area;
33 char* order=NULL;
34
35 /*Domain outline variables: */
36 int nprof;
37 int* profngrids=NULL;
38 double** pprofx=NULL;
39 double** pprofy=NULL;
40 double* xprof=NULL;
41 double* yprof=NULL;
42 int numberofpoints;
43
44 /*Rift outline variables: */
45 int numrifts;
46 int* riftsnumgrids=NULL;
47 double** riftsgridsx=NULL;
48 double** riftsgridsy=NULL;
49
50 /* error handling: */
51 int noerr=1;
52
53 /* Triangle structures: */
54 struct triangulateio in,out;
55 char options[256];
56
57 /* verify correct usage: */
58 if (nlhs==0 && nrhs==0) {
59 /* special case: */
60 TriMeshUsage();
61 return;
62 }
63
64 if (!( (nlhs==5) ||(nrhs==2) || (nrhs==3) || (nrhs==4) )){
65 mexPrintf(" %s format error.\n", __FUNCT__);
66 TriMeshUsage();
67 printf(" ");
68 mexErrMsgTxt(" ");
69 }
70
71 /*Fetch data needed by Triangle: */
72
73 prhs_counter=0;
74 /*First recover the domain outline file name: */
75 if (!mxIsChar(prhs[prhs_counter])){
76 mexPrintf("%s%s\n",__FUNCT__," error message; first argument should be the domain outline file name!");
77 mexErrMsgTxt(" ");
78 }
79 domainname = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
80 mxGetString(prhs[prhs_counter],domainname,mxGetN(prhs[prhs_counter])+1);
81
82 /*Look for optional rifts file name: */
83 prhs_counter++;
84 if (mxIsChar(prhs[prhs_counter])){
85 riftname = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
86 mxGetString(prhs[prhs_counter],riftname,mxGetN(prhs[prhs_counter])+1);
87 prhs_counter++;
88 }
89
90 /*Recover the mesh density desired:*/
91 area=mxGetScalar(prhs[prhs_counter]);
92
93 /*Optionaly, recover desired order: */
94 prhs_counter++;
95 if (mxIsChar(prhs[prhs_counter])){
96 order = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
97 mxGetString(prhs[prhs_counter],order,mxGetN(prhs[prhs_counter])+1);
98 }
99
100 #ifdef _ISSM_DEBUG_
101 printf("Domain name: %s\n",domainname);
102 printf("Rift name: %s\n",riftname);
103 printf("Maximum area: %16.16lf\n",area);
104 printf("Order: %s\n",order);
105 #endif
106
107
108 /*Start reading the domain outline file: */
109 noerr=DomainOutlineRead(&nprof,&profngrids,&pprofx,&pprofy,domainname);
110 if(!noerr){
111 printf("%s%s%s\n",__FUNCT__," error message reading domain outline ",domainname);
112 mexErrMsgTxt(" ");
113 }
114
115 /*Read rifts file if present: */
116 if(riftname){
117 noerr=DomainOutlineRead(&numrifts,&riftsnumgrids,&riftsgridsx,&riftsgridsy,riftname);
118 if(!noerr){
119 printf("%s%s%s\n",__FUNCT__," error message reading rifts outline ",riftname);
120 mexErrMsgTxt(" ");
121 }
122 }
123
124 #ifdef _ISSM_DEBUG_
125 printf("# points in domain outline: %i\n",nprof);
126 #endif
127 #ifdef _DEBUG2_
128 for (i=0;i<nprof;i++){
129 printf("Profile #%i\n",i);
130 x=pprofx[i];
131 y=pprofy[i];
132 for (j=0;j<profngrids[i];j++){
133 printf(" %lf %lf\n",x[j],y[j]);
134 }
135 }
136 #endif
137
138
139 #ifdef _ISSM_DEBUG_
140 printf("# points in rift outline: %i\n",numrifts);
141 #endif
142 #ifdef _DEBUG2_
143 for (i=0;i<numrifts;i++){
144 printf("Profile #%i\n",i);
145 x=riftsgridsx[i];
146 y=riftsgridsy[i];
147 for (j=0;j<riftsnumgrids[i];j++){
148 printf(" %lf %lf\n",x[j],y[j]);
149 }
150 }
151 #endif
152
153 /*Create initial triangulation to call triangulate():*/
154 numberofpoints=0;
155 for (i=0;i<nprof;i++){
156 numberofpoints+=profngrids[i];
157 }
158 if (riftname){
159 for (i=0;i<numrifts;i++){
160 numberofpoints+=riftsnumgrids[i];
161 }
162 }
163 in.numberofpoints=numberofpoints;
164
165 in.numberofpointattributes=1;
166 in.pointlist = (REAL *) mxMalloc(in.numberofpoints * 2 * sizeof(REAL));
167
168 counter=0;
169 for (i=0;i<nprof;i++){
170 xprof=pprofx[i];
171 yprof=pprofy[i];
172 for (j=0;j<profngrids[i];j++){
173 in.pointlist[2*counter+0]=xprof[j];
174 in.pointlist[2*counter+1]=yprof[j];
175 counter++;
176 }
177 }
178 if(riftname){
179 for (i=0;i<numrifts;i++){
180 xprof=riftsgridsx[i];
181 yprof=riftsgridsy[i];
182 for (j=0;j<riftsnumgrids[i];j++){
183 in.pointlist[2*counter+0]=xprof[j];
184 in.pointlist[2*counter+1]=yprof[j];
185 counter++;
186 }
187 }
188 }
189
190 in.pointattributelist = (REAL *) mxMalloc(in.numberofpoints *
191 in.numberofpointattributes *
192 sizeof(REAL));
193 for (i=0;i<in.numberofpoints;i++){
194 in.pointattributelist[i] = 0.0;
195 }
196 in.pointmarkerlist = (int *) mxMalloc(in.numberofpoints * sizeof(int));
197 for(i=0;i<in.numberofpoints;i++){
198 in.pointmarkerlist[i] = 0;
199 }
200
201
202 /*Build segments: */
203 /*Figure out number of segments: holes and closed outlines have as many segments as grids,
204 *for rifts, we have one less segment as we have grids*/
205 in.numberofsegments=0;
206 for (i=0;i<nprof;i++){
207 in.numberofsegments+=profngrids[i];
208 }
209 if (riftname){
210 for (i=0;i<numrifts;i++){
211 in.numberofsegments+=riftsnumgrids[i]-1;
212 }
213 }
214
215 in.segmentlist = (int *) mxMalloc(in.numberofsegments * 2 * sizeof(int));
216 in.segmentmarkerlist = (int *) mxCalloc(in.numberofsegments,sizeof(int));
217 counter=0;
218 backcounter=0;
219 for (i=0;i<nprof;i++){
220 for (j=0;j<(profngrids[i]-1);j++){
221 in.segmentlist[2*counter+0]=counter;
222 in.segmentlist[2*counter+1]=counter+1;
223 in.segmentmarkerlist[counter]=0;
224 counter++;
225 }
226 /*Close this profile: */
227 in.segmentlist[2*counter+0]=counter;
228 in.segmentlist[2*counter+1]=backcounter;
229 in.segmentmarkerlist[counter]=0;
230 counter++;
231 backcounter=counter;
232 }
233 counter2=counter;
234 if(riftname){
235 for (i=0;i<numrifts;i++){
236 for (j=0;j<(riftsnumgrids[i]-1);j++){
237 in.segmentlist[2*counter2+0]=counter;
238 in.segmentlist[2*counter2+1]=counter+1;
239 in.segmentmarkerlist[counter2]=2+i;
240 counter2++;
241 counter++;
242 }
243 counter++;
244 }
245 }
246
247
248 /*Build regions: */
249 in.numberofregions = 0;
250
251 /*Build holes: */
252 in.numberofholes = nprof-1; /*everything is a hole, but for the first profile.*/
253 in.holelist = (REAL *) mxMalloc(in.numberofholes * 2 * sizeof(REAL));
254 for (i=0;i<nprof-1;i++){
255 /*We are looking for a grid that lies inside the hole: */
256 GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],profngrids[i+1],pprofx[i+1],pprofy[i+1]);
257 }
258
259 /* Make necessary initializations so that Triangle can return a */
260 /* triangulation in `out': */
261
262 out.pointlist = (REAL *) NULL;
263 out.pointattributelist = (REAL *) NULL;
264 out.pointmarkerlist = (int *) NULL;
265 out.trianglelist = (int *) NULL;
266 out.triangleattributelist = (REAL *) NULL;
267 out.neighborlist = (int *) NULL;
268 out.segmentlist = (int *) NULL;
269 out.segmentmarkerlist = (int *) NULL;
270 out.edgelist = (int *) NULL;
271 out.edgemarkerlist = (int *) NULL;
272
273 /* Triangulate the points:. Switches are chosen to read and write a */
274 /* PSLG (p), preserve the convex hull (c), number everything from */
275 /* zero (z), assign a regional attribute to each element (A), and */
276 /* produce an edge list (e), a Voronoi diagram (v), and a triangle */
277 /* neighbor list (n). */
278
279 sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
280
281 triangulate(options, &in, &out, NULL);
282 /*report(&out, 0, 1, 1, 1, 1, 0);*/
283
284 /*Allocate index, x and y: */
285 index=(double*)mxMalloc(3*out.numberoftriangles*sizeof(double));
286 x=(double*)mxMalloc(out.numberofpoints*sizeof(double));
287 y=(double*)mxMalloc(out.numberofpoints*sizeof(double));
288 segments=(double*)mxMalloc(3*out.numberofsegments*sizeof(double));
289 segmentmarkerlist=(double*)mxMalloc(out.numberofsegments*sizeof(double));
290
291 for (i = 0; i < out.numberoftriangles; i++) {
292 for (j = 0; j < out.numberofcorners; j++) {
293 *(index+3*i+j)=(double)out.trianglelist[i * out.numberofcorners + j]+1;
294 }
295 }
296 for (i = 0; i < out.numberofpoints; i++) {
297 x[i]=out.pointlist[i * 2 + 0];
298 y[i]=out.pointlist[i * 2 + 1];
299 }
300
301 for (i = 0; i < out.numberofsegments; i++) {
302 segments[3*i+0]=(double)out.segmentlist[i*2+0]+1;
303 segments[3*i+1]=(double)out.segmentlist[i*2+1]+1;
304 segmentmarkerlist[i]=(double)out.segmentmarkerlist[i];
305 }
306
307 /*Associate elements with segments: */
308 noerr=AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
309
310 /*Order segments so that their normals point outside the domain: */
311 if(!strcmp(order,"yes")){
312 OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
313 }
314
315 /*Output : */
316 pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
317 mxSetM(pmxa_array,3);
318 mxSetN(pmxa_array,out.numberoftriangles);
319 mxSetPr(pmxa_array,index);
320 mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "'");
321
322 pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
323 mxSetM(pmxa_array,1);
324 mxSetN(pmxa_array,out.numberofpoints);
325 mxSetPr(pmxa_array,x);
326 mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "'");
327
328 pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
329 mxSetM(pmxa_array,1);
330 mxSetN(pmxa_array,out.numberofpoints);
331 mxSetPr(pmxa_array,y);
332 mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "'");
333
334 pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
335 mxSetM(pmxa_array,3);
336 mxSetN(pmxa_array,out.numberofsegments);
337 mxSetPr(pmxa_array,segments);
338 mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "'");
339
340 pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
341 mxSetM(pmxa_array,1);
342 mxSetN(pmxa_array,out.numberofsegments);
343 mxSetPr(pmxa_array,segmentmarkerlist);
344 mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "'");
345
346 return;
347}
348
349void TriMeshUsage(void)
350{
351 printf("\n");
352 printf(" usage: [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,riftsoutlinename,area,ordered) \n");
353 printf(" where: index,x,y defines a triangulation, segments is an array made \n");
354 printf(" of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment \n");
355 printf(" (if rifts are present, markers >=2 flag them ), outlinefilename an Argus domain outline file.\n");
356 printf(" riftsoutlinename is an Argus domain file, defining rifts (ie: open profiles), \n");
357 printf(" area is the maximum area desired for any element of the resulting mesh. \n");
358 printf(" and ordered is a string ('yes' or 'no') that determines whether segments are output in the \n");
359 printf(" order they are made by Triangle (ie none), or ordered counter clockwise around the domain outline.\n");
360 printf(" riftsoutlinename and ordered are optional arguments.\n");
361 printf("\n");
362}
363
364
365
366
Note: See TracBrowser for help on using the repository browser.