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

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

Fixed some bugs in DakotaResponses and InputUpdateFromDakota modules.
Took out noerr variable in many of the modules.
File size: 9.5 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, "'");
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, "'");
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, "'");
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, "'");
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, "'");
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.