[1] | 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 |
|
---|
| 10 | void 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;
|
---|
[8306] | 37 | int* profnvertices=NULL;
|
---|
[1] | 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;
|
---|
[8306] | 46 | int* riftsnumvertices=NULL;
|
---|
| 47 | double** riftsverticesx=NULL;
|
---|
| 48 | double** riftsverticesy=NULL;
|
---|
[1] | 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: */
|
---|
[8737] | 98 | if(!DomainOutlineRead(&nprof,&profnvertices,&pprofx,&pprofy,NULL,domainname)){
|
---|
[1] | 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){
|
---|
[8737] | 105 | if(!DomainOutlineRead(&numrifts,&riftsnumvertices,&riftsverticesx,&riftsverticesy,NULL,riftname)){
|
---|
[1] | 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++){
|
---|
[8306] | 114 | numberofpoints+=profnvertices[i];
|
---|
[1] | 115 | }
|
---|
| 116 | if (riftname){
|
---|
| 117 | for (i=0;i<numrifts;i++){
|
---|
[8306] | 118 | numberofpoints+=riftsnumvertices[i];
|
---|
[1] | 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];
|
---|
[8306] | 130 | for (j=0;j<profnvertices[i];j++){
|
---|
[1] | 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++){
|
---|
[8306] | 138 | xprof=riftsverticesx[i];
|
---|
| 139 | yprof=riftsverticesy[i];
|
---|
| 140 | for (j=0;j<riftsnumvertices[i];j++){
|
---|
[1] | 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: */
|
---|
[8306] | 160 | /*Figure out number of segments: holes and closed outlines have as many segments as vertices,
|
---|
| 161 | *for rifts, we have one less segment as we have vertices*/
|
---|
[1] | 162 | in.numberofsegments=0;
|
---|
| 163 | for (i=0;i<nprof;i++){
|
---|
[8306] | 164 | in.numberofsegments+=profnvertices[i];
|
---|
[1] | 165 | }
|
---|
| 166 | if (riftname){
|
---|
| 167 | for (i=0;i<numrifts;i++){
|
---|
[8306] | 168 | in.numberofsegments+=riftsnumvertices[i]-1;
|
---|
[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++){
|
---|
[8306] | 177 | for (j=0;j<(profnvertices[i]-1);j++){
|
---|
[1] | 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++){
|
---|
[8306] | 193 | for (j=0;j<(riftsnumvertices[i]-1);j++){
|
---|
[1] | 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++){
|
---|
[8306] | 212 | /*We are looking for a vertex that lies inside the hole: */
|
---|
| 213 | GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],profnvertices[i+1],pprofx[i+1],pprofy[i+1]);
|
---|
[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 |
|
---|
[1130] | 236 | sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
|
---|
[1] | 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: */
|
---|
[4453] | 265 | AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
|
---|
[1] | 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);
|
---|
[4526] | 277 | mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "transpose");
|
---|
[1] | 278 |
|
---|
| 279 | pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
|
---|
| 280 | mxSetM(pmxa_array,1);
|
---|
| 281 | mxSetN(pmxa_array,out.numberofpoints);
|
---|
| 282 | mxSetPr(pmxa_array,x);
|
---|
[4526] | 283 | mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "transpose");
|
---|
[1] | 284 |
|
---|
| 285 | pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
|
---|
| 286 | mxSetM(pmxa_array,1);
|
---|
| 287 | mxSetN(pmxa_array,out.numberofpoints);
|
---|
| 288 | mxSetPr(pmxa_array,y);
|
---|
[4526] | 289 | mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "transpose");
|
---|
[1] | 290 |
|
---|
| 291 | pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
|
---|
| 292 | mxSetM(pmxa_array,3);
|
---|
| 293 | mxSetN(pmxa_array,out.numberofsegments);
|
---|
| 294 | mxSetPr(pmxa_array,segments);
|
---|
[4526] | 295 | mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "transpose");
|
---|
[1] | 296 |
|
---|
| 297 | pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
|
---|
| 298 | mxSetM(pmxa_array,1);
|
---|
| 299 | mxSetN(pmxa_array,out.numberofsegments);
|
---|
| 300 | mxSetPr(pmxa_array,segmentmarkerlist);
|
---|
[4526] | 301 | mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "transpose");
|
---|
[1] | 302 |
|
---|
| 303 | return;
|
---|
| 304 | }
|
---|
| 305 |
|
---|
| 306 | void 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 |
|
---|