Changeset 11882

Timestamp:

04/03/12 16:20:08 (13 years ago)

Author:

Mathieu Morlighem

Message:

Back to normal for now

File:

• issm/trunk-jpl/src/mex/TriMesh/TriMesh.cpp (modified) (2 diffs)

issm/trunk-jpl/src/mex/TriMesh/TriMesh.cpp

@@ -1 +1 @@
 /*
- * TriMesh: mesh a domain using an .exp file
+ * TriMesh: out of a domain outline file ( Argus format ), 
+ * use the Triangle package to create a triangular mesh 
+ *
  */
 
 #include "./TriMesh.h"
 
-void mexFunction(       int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[] ){
+
+void mexFunction(       int nlhs, mxArray* plhs[],
+                                        int nrhs, const mxArray* prhs[] )
+{
+
+
+        /*Matlab arrays: */
+        mxArray* pmxa_array=NULL;
+        int i,j;
+        int counter,counter2,backcounter;
+        int prhs_counter;
+        
+        /* returned quantities: */
+
+        double* index=NULL;
+        double* x=NULL;
+        double* y=NULL;
+        double* segments=NULL;
+        double*    segmentmarkerlist=NULL;
 
         /* input: */
         char*  domainname=NULL;
+        char*  riftname=NULL;
         double area;
-        bool   order;
-
-        /*intermediary: */
-        DataSet* domain=NULL;
-
-        /* output: */
-        Matrix* index=NULL;
-        Vector* x=NULL;
-        Vector* y=NULL;
-        Matrix* segments=NULL;
-        Vector* segmentmarkerlist=NULL;
-
-        /*Boot module: */
-        MODULEBOOT();
-printf("ok1\n");
-
-        /*checks on arguments on the matlab side: */
-        CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&TriMeshUsage);
-printf("ok2\n");
-
-        /*Fetch data needed for meshing: */
-
-        FetchMatlabData(&domainname,DOMAINOUTLINE);
-        FetchMatlabData(&area,AREA);
-        FetchMatlabData(&order,ORDER);
-printf("ok3 %g %s\n",area,domainname);
-
-        /*Read domain outline: */
-        domain=DomainOutlineRead(domainname,false);
-printf("ok4\n");
-
-        /*call x core: */
-        TriMeshx(&index,&x,&y,&segments,&segmentmarkerlist,domain,area,order);
-printf("ok5\n");
-
-        /*write outputs: */
-        /*WriteMatlabData(INDEX,index);
-        WriteMatlabData(X,x);
-        WriteMatlabData(Y,y);
-        WriteMatlabData(SEGMENTS,segments);
-        WriteMatlabData(SEGMENTMARKERLIST,segmentmarkerlist);*/
-
-        /*free ressources: */
-        delete domain;
-        /*
-        xdelete(&index);
-        xdelete(&x);
-        xdelete(&y);
-        xdelete(&segments);
-        xdelete(&segmentmarkerlist);*/
-
-        /*end module: */
-        MODULEEND();
+        char*  order=NULL;
+        
+        /*Domain outline variables: */
+        int      nprof;
+        int*     profnvertices=NULL;
+        double** pprofx=NULL;
+        double** pprofy=NULL;
+        double*  xprof=NULL;
+        double*  yprof=NULL;
+        int      numberofpoints;
+
+        /*Rift outline variables: */
+        int      numrifts;
+        int*     riftsnumvertices=NULL;
+        double** riftsverticesx=NULL;
+        double** riftsverticesy=NULL;
+
+        /* Triangle structures: */
+        struct triangulateio in,out;
+        char   options[256];
+
+        /* verify correct usage: */
+        if (nlhs==0 && nrhs==0) {
+                /* special case: */
+                TriMeshUsage();
+                return;
+        }
+
+        if (!(  (nlhs==5) ||(nrhs==2) || (nrhs==3)  || (nrhs==4) )){
+                mexPrintf("   %s format error.\n", __FUNCT__);
+                TriMeshUsage();
+                printf("   ");
+                mexErrMsgTxt(" ");
+        }
+
+        /*Fetch data needed by Triangle: */
+
+        prhs_counter=0;
+        /*First recover the domain outline file name: */
+        if (!mxIsChar(prhs[prhs_counter])){
+                mexPrintf("%s%s\n",__FUNCT__," error message; first argument should be the domain outline file name!");
+                mexErrMsgTxt(" ");
+        }
+        domainname = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
+        mxGetString(prhs[prhs_counter],domainname,mxGetN(prhs[prhs_counter])+1);
+
+        /*Look for optional rifts file name: */
+        prhs_counter++;
+        if (mxIsChar(prhs[prhs_counter])){
+                riftname = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
+                mxGetString(prhs[prhs_counter],riftname,mxGetN(prhs[prhs_counter])+1);
+                prhs_counter++;
+        }
+
+        /*Recover the mesh density desired:*/
+        area=mxGetScalar(prhs[prhs_counter]);
+
+        /*Optionaly, recover desired order: */
+        prhs_counter++;
+        if (mxIsChar(prhs[prhs_counter])){
+                order = (char *) mxMalloc((mxGetN(prhs[prhs_counter])+1)*sizeof(char));
+                mxGetString(prhs[prhs_counter],order,mxGetN(prhs[prhs_counter])+1);
+        }
+        
+        /*Start reading the domain outline file: */
+        if(!DomainOutlineRead(&nprof,&profnvertices,&pprofx,&pprofy,NULL,domainname,false)){
+                printf("%s%s%s\n",__FUNCT__," error message reading domain outline ",domainname);
+                mexErrMsgTxt(" ");
+        }
+
+        /*Read rifts file if present: */
+        if(riftname){
+                if(!DomainOutlineRead(&numrifts,&riftsnumvertices,&riftsverticesx,&riftsverticesy,NULL,riftname,false)){
+                        printf("%s%s%s\n",__FUNCT__," error message reading rifts outline ",riftname);
+                        mexErrMsgTxt(" ");
+                }
+        }
+
+        /*Create initial triangulation to call triangulate():*/
+        numberofpoints=0;
+        for (i=0;i<nprof;i++){
+                numberofpoints+=profnvertices[i];
+        }
+        if (riftname){
+                for (i=0;i<numrifts;i++){
+                        numberofpoints+=riftsnumvertices[i];
+                }
+        }
+        in.numberofpoints=numberofpoints;
+
+        in.numberofpointattributes=1;
+        in.pointlist = (REAL *) mxMalloc(in.numberofpoints * 2 * sizeof(REAL));
+
+        counter=0;
+        for (i=0;i<nprof;i++){
+                xprof=pprofx[i];
+                yprof=pprofy[i];
+                for (j=0;j<profnvertices[i];j++){
+                        in.pointlist[2*counter+0]=xprof[j];
+                        in.pointlist[2*counter+1]=yprof[j];
+                        counter++;
+                }
+        }
+        if(riftname){
+                for (i=0;i<numrifts;i++){
+                        xprof=riftsverticesx[i];
+                        yprof=riftsverticesy[i];
+                        for (j=0;j<riftsnumvertices[i];j++){
+                                in.pointlist[2*counter+0]=xprof[j];
+                                in.pointlist[2*counter+1]=yprof[j];
+                                counter++;
+                        }
+                }
+        }
+        
+        in.pointattributelist = (REAL *) mxMalloc(in.numberofpoints *
+                                                                                  in.numberofpointattributes *
+                                                                                  sizeof(REAL));
+        for (i=0;i<in.numberofpoints;i++){
+                in.pointattributelist[i] = 0.0;
+        }
+        in.pointmarkerlist = (int *) mxMalloc(in.numberofpoints * sizeof(int));
+        for(i=0;i<in.numberofpoints;i++){
+                in.pointmarkerlist[i] = 0;
+        }
+
+        /*Build segments: */
+        /*Figure out number of segments: holes and closed outlines have as many segments as vertices, 
+         *for rifts, we have one less segment as we have vertices*/
+        in.numberofsegments=0;
+        for (i=0;i<nprof;i++){
+                in.numberofsegments+=profnvertices[i];
+        }
+        if (riftname){
+                for (i=0;i<numrifts;i++){
+                        in.numberofsegments+=riftsnumvertices[i]-1;
+                }
+        }
+        
+        in.segmentlist = (int *) mxMalloc(in.numberofsegments * 2 * sizeof(int));
+        in.segmentmarkerlist = (int *) mxCalloc(in.numberofsegments,sizeof(int));
+        counter=0;
+        backcounter=0;
+        for (i=0;i<nprof;i++){
+                for (j=0;j<(profnvertices[i]-1);j++){
+                        in.segmentlist[2*counter+0]=counter;
+                        in.segmentlist[2*counter+1]=counter+1;
+                        in.segmentmarkerlist[counter]=0;
+                        counter++;
+                }
+                /*Close this profile: */
+                 in.segmentlist[2*counter+0]=counter;
+                 in.segmentlist[2*counter+1]=backcounter;
+                 in.segmentmarkerlist[counter]=0;
+                 counter++;
+                 backcounter=counter;
+        }
+        counter2=counter;
+        if(riftname){
+                for (i=0;i<numrifts;i++){
+                        for (j=0;j<(riftsnumvertices[i]-1);j++){
+                                in.segmentlist[2*counter2+0]=counter;
+                                in.segmentlist[2*counter2+1]=counter+1;
+                                in.segmentmarkerlist[counter2]=2+i;
+                                counter2++;
+                                counter++;
+                        }
+                        counter++;
+                }
+        }
+
+        
+        /*Build regions: */
+        in.numberofregions = 0;
+
+        /*Build holes: */
+        in.numberofholes = nprof-1; /*everything is a hole, but for the first profile.*/
+        in.holelist = (REAL *) mxMalloc(in.numberofholes * 2 * sizeof(REAL));
+        for (i=0;i<nprof-1;i++){
+                /*We are looking for a vertex that lies inside the hole: */
+                GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],profnvertices[i+1],pprofx[i+1],pprofy[i+1]);
+        }
+
+        /* Make necessary initializations so that Triangle can return a */
+        /*   triangulation in `out': */
+
+        out.pointlist = (REAL *) NULL;            
+        out.pointattributelist = (REAL *) NULL;
+        out.pointmarkerlist = (int *) NULL; 
+        out.trianglelist = (int *) NULL;          
+        out.triangleattributelist = (REAL *) NULL;
+        out.neighborlist = (int *) NULL;         
+        out.segmentlist = (int *) NULL;
+        out.segmentmarkerlist = (int *) NULL;
+        out.edgelist = (int *) NULL;             
+        out.edgemarkerlist = (int *) NULL;   
+
+        /* Triangulate the points:.  Switches are chosen to read and write a  */
+        /*   PSLG (p), preserve the convex hull (c), number everything from  */
+        /*   zero (z), assign a regional attribute to each element (A), and  */
+        /*   produce an edge list (e), a Voronoi diagram (v), and a triangle */
+        /*   neighbor list (n).                                              */
+
+        sprintf(options,"%s%lf","pQzDq30ia",area); /*replace V by Q to quiet down the logging*/
+  
+        triangulate(options, &in, &out, NULL);
+        /*report(&out, 0, 1, 1, 1, 1, 0);*/
+
+        /*Allocate index, x and y: */
+        index=(double*)mxMalloc(3*out.numberoftriangles*sizeof(double));
+        x=(double*)mxMalloc(out.numberofpoints*sizeof(double));
+        y=(double*)mxMalloc(out.numberofpoints*sizeof(double));
+        segments=(double*)mxMalloc(3*out.numberofsegments*sizeof(double));
+        segmentmarkerlist=(double*)mxMalloc(out.numberofsegments*sizeof(double));
+
+        for (i = 0; i < out.numberoftriangles; i++) {
+                for (j = 0; j < out.numberofcorners; j++) {
+                        *(index+3*i+j)=(double)out.trianglelist[i * out.numberofcorners + j]+1;
+                }
+        }
+        for (i = 0; i < out.numberofpoints; i++) {
+                x[i]=out.pointlist[i * 2 + 0];
+                y[i]=out.pointlist[i * 2 + 1];
+        }
+        
+        for (i = 0; i < out.numberofsegments; i++) {
+                segments[3*i+0]=(double)out.segmentlist[i*2+0]+1;
+                segments[3*i+1]=(double)out.segmentlist[i*2+1]+1;
+                segmentmarkerlist[i]=(double)out.segmentmarkerlist[i];
+        }
+
+        /*Associate elements with segments: */
+        AssociateSegmentToElement(&segments,out.numberofsegments,index,out.numberoftriangles);
+
+        /*Order segments so that their normals point outside the domain: */
+        if(!strcmp(order,"yes")){
+                OrderSegments(&segments,out.numberofsegments, index,out.numberoftriangles);
+        }
+
+        /*Output : */
+        pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+        mxSetM(pmxa_array,3);
+        mxSetN(pmxa_array,out.numberoftriangles);
+        mxSetPr(pmxa_array,index);
+        mexCallMATLAB( 1, &plhs[0], 1, &pmxa_array, "transpose");
+
+        pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+        mxSetM(pmxa_array,1);
+        mxSetN(pmxa_array,out.numberofpoints);
+        mxSetPr(pmxa_array,x);
+        mexCallMATLAB( 1, &plhs[1], 1, &pmxa_array, "transpose");
+
+        pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+        mxSetM(pmxa_array,1);
+        mxSetN(pmxa_array,out.numberofpoints);
+        mxSetPr(pmxa_array,y);
+        mexCallMATLAB( 1, &plhs[2], 1, &pmxa_array, "transpose");
+
+        pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+        mxSetM(pmxa_array,3);
+        mxSetN(pmxa_array,out.numberofsegments);
+        mxSetPr(pmxa_array,segments);
+        mexCallMATLAB( 1, &plhs[3], 1, &pmxa_array, "transpose");
+
+        pmxa_array= mxCreateDoubleMatrix(0,0,mxREAL);
+        mxSetM(pmxa_array,1);
+        mxSetN(pmxa_array,out.numberofsegments);
+        mxSetPr(pmxa_array,segmentmarkerlist);
+        mexCallMATLAB( 1, &plhs[4], 1, &pmxa_array, "transpose");
+        
+        return;
 }
 
@@ -68 +307 @@
 {
         printf("\n");
-        printf("   usage: [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,area,ordered) \n");
+        printf("   usage: [index,x,y,segments,segmentmarkers]=TriMesh(domainoutlinefilename,riftsoutlinename,area,ordered) \n");
         printf("      where: index,x,y defines a triangulation, segments is an array made \n");
-        printf("      of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment, \n");
-        printf("      outlinefilename an Argus domain outline file, \n");
-        printf("      area is the maximum area desired for any element of the resulting mesh, \n");
-        printf("      and ordered is a bool that determines whether segments are output in the \n");
+        printf("      of exterior segments to the mesh domain outline, segmentmarkers is an array flagging each segment \n");
+        printf("      (if rifts are present, markers >=2 flag them ), outlinefilename an Argus domain outline file.\n");
+        printf("      riftsoutlinename is an Argus domain file, defining rifts (ie: open profiles), \n");
+        printf("      area is the maximum area desired for any element of the resulting mesh. \n");
+        printf("      and ordered is a string ('yes' or 'no') that determines whether segments are output in the \n");
         printf("      order they are made by Triangle (ie none), or ordered counter clockwise around the domain outline.\n");
+        printf("      riftsoutlinename and ordered are optional arguments.\n");
         printf("\n");
 }
+
+
+
+