Changeset 12071

Timestamp:

04/20/12 08:39:29 (13 years ago)

Author:

Mathieu Morlighem

Message:

Improved TriMeshRifts, to be merged to TriMesh

Location:

issm/trunk-jpl/src/modules

Files:

• TriMeshProcessRifts/TriMeshProcessRifts.cpp (modified) (1 diff)
• TriMeshRifts/TriMeshRifts.cpp (modified) (1 diff)
• TriMeshRifts/TriMeshRifts.h (modified) (1 diff)

issm/trunk-jpl/src/modules/TriMeshProcessRifts/TriMeshProcessRifts.cpp ¶

@@ -189 +189 @@
                 
                 /*Order rifts so that they start from one tip, go to the other tip, and back: */
-                OrderRifts(&out_riftstips, out_riftssegments,out_riftspairs,numrifts,out_riftsnumsegments,x_in,y_in);
+                OrderRifts(&out_riftstips, out_riftssegments,out_riftspairs,numrifts,out_riftsnumsegments,x_in,y_in,nods,nel);
 
                 /*Create penalty pairs, used by Imp: */

issm/trunk-jpl/src/modules/TriMeshRifts/TriMeshRifts.cpp ¶

@@ -1 +1 @@
 /*
- * TriMeshRifts: out of a domain outline file ( Argus format ), 
- * use the Triangle package to create a triangular mesh 
- *
+ * TriMeshRifts: mesh a domain using an .exp file
  */
 
 #include "./TriMeshRifts.h"
 
-void mexFunction(       int nlhs, mxArray* plhs[],
-                                        int nrhs, const mxArray* prhs[] )
-{
+WRAPPER(TriMeshRifts){
+        
+        /* input: */
+        char   *domainname = NULL;
+        char   *riftsname  = NULL;
+        double  area;
+        bool    order;
 
+        /*intermediary: */
+        DataSet *domain = NULL;
+        DataSet *rifts  = NULL;
 
-        /*Matlab arrays: */
-        mxArray* pmxa_array=NULL;
-        int i,j;
-        int counter,counter2,backcounter;
-        int prhs_counter;
+        /* output: */
+        Matrix *index             = NULL;
+        Vector *x                 = NULL;
+        Vector *y                 = NULL;
+        Matrix *segments          = NULL;
+        Vector *segmentmarkerlist = NULL;
+
+        /*Boot module: */
+        MODULEBOOT();
+
+        /*checks on arguments on the matlab side: */
+        CHECKARGUMENTS(NLHS,NRHS,&TriMeshRiftsUsage);
         
-        /* returned quantities: */
+        /*Fetch data needed for meshing: */
+        FetchData(&domainname,DOMAINOUTLINE);
+        FetchData(&riftsname,RIFTSOUTLINE);
+        FetchData(&area,AREA);
+        FetchData(&order,ORDER);
 
-        double* index=NULL;
-        double* x=NULL;
-        double* y=NULL;
-        double* segments=NULL;
-        double*    segmentmarkerlist=NULL;
+        /*Read domain outline: */
+        domain = DomainOutlineRead(domainname,false);
+        rifts  = DomainOutlineRead(riftsname,false);
 
-        /* input: */
-        char*  domainname=NULL;
-        char*  riftname=NULL;
-        double area;
-        char*  order=NULL;
+        /*call x core: */
+        TriMeshRiftsx(&index,&x,&y,&segments,&segmentmarkerlist,domain,rifts,area,order);
         
-        /*Domain outline variables: */
-        int      nprof;
-        int*     profnvertices=NULL;
-        double** pprofx=NULL;
-        double** pprofy=NULL;
-        double*  xprof=NULL;
-        double*  yprof=NULL;
-        int      numberofpoints;
+        /*write outputs: */
+        WriteData(INDEX,index);
+        WriteData(X,x);
+        WriteData(Y,y);
+        WriteData(SEGMENTS,segments);
+        WriteData(SEGMENTMARKERLIST,segmentmarkerlist);
 
-        /*Rift outline variables: */
-        int      numrifts;
-        int*     riftsnumvertices=NULL;
-        double** riftsverticesx=NULL;
-        double** riftsverticesy=NULL;
+        /*free ressources: */
+        delete domain;
+        delete rifts;
+        xdelete_module(&index);
+        xdelete_module(&x);
+        xdelete_module(&y);
+        xdelete_module(&segments);
+        xdelete_module(&segmentmarkerlist);
 
-        /* Triangle structures: */
-        struct triangulateio in,out;
-        char   options[256];
-
-        /* verify correct usage: */
-        if (nlhs==0 && nrhs==0) {
-                /* special case: */
-                TriMeshRiftsUsage();
-                return;
-        }
-
-        if (!(  (nlhs==5) ||(nrhs==2) || (nrhs==3)  || (nrhs==4) )){
-                mexPrintf("   %s format error.\n", __FUNCT__);
-                TriMeshRiftsUsage();
-                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;
+        /*end module: */
+        MODULEEND();
 }
 
-void TriMeshRiftsUsage(void)
+void TriMeshRiftsUsage(void) //{{{1
 {
         printf("\n");
-        printf("   usage: [index,x,y,segments,segmentmarkers]=TriMeshRifts(domainoutlinefilename,riftsoutlinename,area,ordered) \n");
+        printf("   usage: [index,x,y,segments,segmentmarkers]=TriMeshRifts(domainoutlinefilename,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("      (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("      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("      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");
 }
-
-
-
-
+//}}}

issm/trunk-jpl/src/modules/TriMeshRifts/TriMeshRifts.h ¶

@@ -6 +6 @@
 #define _TRIMESHRIFTS_H_
 
-#include "mex.h"
-#include "triangle.h"
-#include "string.h"
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
 
+/*Very important definition in case we are compiling a python module!: needs to come before header files inclusion*/
+#ifdef _HAVE_PYTHON_
+#define PY_ARRAY_UNIQUE_SYMBOL PythonIOSymbol
+#endif
+
+/*Header files: */
 #include "../../c/include/globals.h"
+#include "../../c/toolkits/toolkits.h"
+#include "../../c/include/include.h"
 #include "../../c/modules/modules.h"
 #include "../../c/Container/Container.h"
 #include "../../c/shared/shared.h"
 #include "../../c/issm-binding.h"
+#include "../../c/io/io.h"
+#include "../../c/EnumDefinitions/EnumDefinitions.h"
 
-void TriMeshRiftsUsage(void);
+#ifdef _HAVE_MATLAB_MODULES_
+/* serial input macros: */
+#define DOMAINOUTLINE  (mxArray *)prhs[0]
+#define RIFTSOUTLINE   (mxArray *)prhs[1]
+#define AREA           (mxArray *)prhs[2]
+#define ORDER          (mxArray *)prhs[3]
+/* serial output macros: */
+#define INDEX             (mxArray**)&plhs[0]
+#define X                 (mxArray**)&plhs[1]
+#define Y                 (mxArray**)&plhs[2]
+#define SEGMENTS          (mxArray**)&plhs[3]
+#define SEGMENTMARKERLIST (mxArray**)&plhs[4]
+#endif
+
+#ifdef _HAVE_PYTHON_MODULES_
+/* serial input macros: */
+#define DOMAINOUTLINE PyTuple_GetItem(args,0)
+#define RIFTSOUTLINE  PyTuple_GetItem(args,1)
+#define AREA          PyTuple_GetItem(args,2)
+#define ORDER         PyTuple_GetItem(args,3)
+/* serial output macros: */
+#define INDEX             output,0
+#define X                 output,1
+#define Y                 output,2
+#define SEGMENTS          output,3
+#define SEGMENTMARKERLIST output,4
+#endif
 
 #undef __FUNCT__ 
 #define __FUNCT__ "TriMeshRifts"
 
-#endif
+/* serial arg counts: */
+#undef NLHS
+#define NLHS  5
+#undef NRHS
+#define NRHS  4
+
+/* local prototypes: */
+void TriMeshRiftsUsage(void);
+
+#endif  /* _TRIMESHRIFTS_H_*/