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

Last change on this file since 4526 was 4526, checked in by Mathieu Morlighem, 15 years ago

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