Changeset 5032

Timestamp:

08/06/10 14:48:41 (15 years ago)

Author:

Eric.Larour

Message:

Added InterpFromMesh2d module, from ISSM2.0, to go around some of Bamg
error messages. We'll get rid of it later, when bugs are addressed.

Location:

issm/trunk/src

Files:

• c/Makefile.am (modified) (2 diffs)
• c/modules/InterpFromMesh2dx (added)
• c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp (added)
• c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h (added)
• c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (modified) (5 diffs)
• c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h (modified) (1 diff)
• c/modules/modules.h (modified) (1 diff)
• mex/InterpFromMesh2d (added)
• mex/InterpFromMesh2d/1 (added)
• mex/InterpFromMesh2d/InterpFromMesh2d.cpp (added)
• mex/InterpFromMesh2d/InterpFromMesh2d.h (added)
• mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp (modified) (7 diffs)
• mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.h (modified) (2 diffs)
• mex/Makefile.am (modified) (2 diffs)

issm/trunk/src/c/Makefile.am ¶

@@ -479 +479 @@
                                         ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
                                         ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
+                                        ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\
+                                        ./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\
                                         ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
                                         ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h\
@@ -1015 +1017 @@
                                         ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
                                         ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
+                                        ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\
+                                        ./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\
                                         ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
                                         ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\

issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp ¶

@@ -8 +8 @@
 #include "../../toolkits/toolkits.h"
 #include "../../objects/objects.h"
+#include "../modules.h"
 
 using namespace bamg;
@@ -13 +14 @@
 
 int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
-                        double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value){
-
+                        double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double* default_values,int num_default_values, Contour** contours, int numcontours){
+        
         /*Output*/
         double* data_interp=NULL;
@@ -30 +31 @@
         int verbose=0;
 
+        /*default values: */
+        Vec    vec_incontour=NULL;
+        double*    incontour=NULL;
+        bool   skip_bamg=false;
+
+
         /*Checks*/
         if (data_cols<=0){
@@ -36 +43 @@
         if (data_rows!=nods_data && data_rows!=nels_data){
                 ISSMERROR("data provided should have either %i or %i lines (not %i)",nods_data,nels_data,data_rows);
+        }
+        if((num_default_values) && (data_cols>1)){
+                ISSMERROR("data provided can only have 1 column if a default value is provided");
+        }
+        
+        /*If default values supplied, figure out which nodes are inside the contour, including the border of the contour: */
+        if(num_default_values){
+                ContourToNodesx( &vec_incontour,x_interp,y_interp,nods_interp,contours,numcontours,1);
+                VecToMPISerial(&incontour,vec_incontour);
         }
 
@@ -50 +66 @@
         if (verbose) printf("Loop over the nodes\n");
         for(i=0;i<nods_interp;i++){
+                
+                /*reset skip_bamg: */
+                skip_bamg=false;
 
-                //Get current point coordinates
-                r.x=x_interp[i]; r.y=y_interp[i];
-                I2 I=Th.toI2(r);
-
-                //Find triangle holding r/I
-                Triangle &tb=*Th.FindTriangleContaining(I,dete);
-
-                // internal point 
-                if (tb.det>0){ 
-                        //Area coordinate
-                        areacoord[0]= (double) dete[0]/ tb.det;
-                        areacoord[1]= (double) dete[1] / tb.det;
-                        areacoord[2]= (double) dete[2] / tb.det;
-                        //3 vertices of the triangle
-                        i0=Th.Number(tb[0]);
-                        i1=Th.Number(tb[1]);
-                        i2=Th.Number(tb[2]);
-                        //triangle number
-                        it=Th.Number(tb);
+                /*figure out if we should skip bamg logic: */
+                if(num_default_values){
+                        if(!incontour[i]){
+                                /*This node is not inside the contour. Skip Bamg logic and apply default value.: */
+                                skip_bamg=true;
+                        }
                 }
 
-                //external point
-                else {
-                        //Get closest adjacent triangle (inside the mesh)
-                        TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
-                        int k=ta;
-                        Triangle &tc=*(Triangle*)ta;
-                        //Area coordinate
-                        areacoord[VerticesOfTriangularEdge[k][1]] = aa;
-                        areacoord[VerticesOfTriangularEdge[k][0]] = bb;
-                        areacoord[OppositeVertex[k]] = 1 - aa -bb;
-                        //3 vertices of the triangle
-                        i0=Th.Number(tc[0]);
-                        i1=Th.Number(tc[1]);
-                        i2=Th.Number(tc[2]);
-                        //triangle number
-                        it=Th.Number(tc);
-                }
+                if(skip_bamg==false){
 
-                /*Last step, P1 interpolation*/
-                if (data_rows==nods_data){
-                        for (j=0;j<data_cols;j++){
-                                data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j];
+                        //Get current point coordinates
+                        r.x=x_interp[i]; r.y=y_interp[i];
+                        I2 I=Th.toI2(r);
+
+                        //Find triangle holding r/I
+                        Triangle &tb=*Th.FindTriangleContaining(I,dete);
+
+                        // internal point 
+                        if (tb.det>0){ 
+                                //Area coordinate
+                                areacoord[0]= (double) dete[0]/ tb.det;
+                                areacoord[1]= (double) dete[1] / tb.det;
+                                areacoord[2]= (double) dete[2] / tb.det;
+                                //3 vertices of the triangle
+                                i0=Th.Number(tb[0]);
+                                i1=Th.Number(tb[1]);
+                                i2=Th.Number(tb[2]);
+                                //triangle number
+                                it=Th.Number(tb);
+                        }
+                        //external point
+                        else {
+                                //Get closest adjacent triangle (inside the mesh)
+                                TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
+                                int k=ta;
+                                Triangle &tc=*(Triangle*)ta;
+                                //Area coordinate
+                                areacoord[VerticesOfTriangularEdge[k][1]] = aa;
+                                areacoord[VerticesOfTriangularEdge[k][0]] = bb;
+                                areacoord[OppositeVertex[k]] = 1 - aa -bb;
+                                //3 vertices of the triangle
+                                i0=Th.Number(tc[0]);
+                                i1=Th.Number(tc[1]);
+                                i2=Th.Number(tc[2]);
+                                //triangle number
+                                it=Th.Number(tc);
+                        }
+                        
+                        if (data_rows==nods_data){
+                                for (j=0;j<data_cols;j++){
+                                        data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j];
+                                }
+                        }
+                        else{
+                                for (j=0;j<data_cols;j++){
+                                        if (it<0 || it>=nels_data){
+                                                ISSMERROR("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data);
+                                        }
+                                        data_interp[i*data_cols+j]=data[data_cols*it+j];
+                                }
                         }
                 }
                 else{
-                        for (j=0;j<data_cols;j++){
-                                if (it<0 || it>=nels_data){
-                                        ISSMERROR("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data);
-                                }
-                                data_interp[i*data_cols+j]=data[data_cols*it+j];
-                        }
+                        if(num_default_values==1) data_interp[i]=default_values[0];
+                        else data_interp[i]=default_values[i];
                 }
-
         }
 

issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h ¶

@@ -10 +10 @@
 /* local prototypes: */
 int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
-                        double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value);
+                        double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double* default_values,int num_default_values,Contour** contours, int numcontours );
 
 #endif

issm/trunk/src/c/modules/modules.h ¶

@@ -30 +30 @@
 #include "./InputDuplicatex/InputDuplicatex.h"
 #include "./InputScalex/InputScalex.h"
+#include "./InterpFromMesh2dx/InterpFromMesh2dx.h"
 #include "./InterpFromGridToMeshx/InterpFromGridToMeshx.h"
 #include "./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h"

issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp ¶

@@ -5 +5 @@
 
 void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
+
+        int i;
 
         /*input: */
@@ -26 +28 @@
         int     y_interp_rows;
 
-        double default_value;
+        double* default_values=NULL;
+        int     num_default_values=0;
+
+        //contours
+        mxArray*  matlabstructure=NULL;
+        Contour** contours=NULL;
+        int       numcontours;
+        Contour*  contouri=NULL;
 
         /*Intermediary*/
@@ -41 +50 @@
 
         /*checks on arguments on the matlab side: */
-        CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InterpFromMeshToMesh2dUsage);
+        if(nlhs!=NLHS){
+                InterpFromMeshToMesh2dUsage();
+                ISSMERROR("InterpFromMeshToMesh2dUsage usage error");
+        }
+        if((nrhs!=6) & (nrhs!=8)){
+                InterpFromMeshToMesh2dUsage();
+                ISSMERROR("InterpFromMeshToMesh2dUsage usage error");
+        }
 
         /*Input datasets: */
@@ -51 +67 @@
         FetchData(&x_interp,&x_interp_rows,NULL,XINTERPHANDLE);
         FetchData(&y_interp,&y_interp_rows,NULL,YINTERPHANDLE);
-        FetchData(&default_value,DEFAULTHANDLE);
+
+        if(nrhs==8){
+                
+                /*Call expread on filename to build a contour array in the matlab workspace: */
+                mexCallMATLAB( 1, &matlabstructure, 1, (mxArray**)&FILENAME, "expread");
+
+                /*default values: */
+                FetchData(&default_values,&num_default_values,DEFAULTHANDLE);
+
+                /*contours: */
+                numcontours=mxGetNumberOfElements(matlabstructure);
+                contours=(Contour**)xmalloc(numcontours*sizeof(Contour*));
+                for(i=0;i<numcontours;i++){
+                        //allocate
+                        contouri=(Contour*)xmalloc(sizeof(Contour));
+                        //retrieve dimension of this contour.
+                        contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));
+                        //set pointers.
+                        contouri->x=mxGetPr(mxGetField(matlabstructure,i,"x"));
+                        contouri->y=mxGetPr(mxGetField(matlabstructure,i,"y"));
+                        *(contours+i)=contouri;
+                }
+
+                /* Debugging of contours :{{{1*/
+                /*for(i=0;i<numcontours;i++){
+                  printf("\nContour echo: contour number  %i / %i\n",i+1,numcontours);
+                  contouri=*(contours+i);
+                  printf("   Number of grids %i\n",contouri->nods);
+                  for (j=0;j<contouri->nods;j++){
+                  printf("   %lf %lf\n",*(contouri->x+j),*(contouri->y+j));
+                  }
+                  }*/
+                /*}}}*/
+        }
+        else{
+                num_default_values=0;
+                default_values=NULL;
+                numcontours=0;
+                contours=NULL;
+        }
 
         /*some checks*/
@@ -71 +126 @@
         /* Run core computations: */
         if (verbose) printf("Call core\n");
-        InterpFromMeshToMesh2dx(&data_interp,index,x_data,y_data,nods_data,nels_data,data,data_rows,data_cols,x_interp,y_interp,nods_interp,default_value);
+        InterpFromMeshToMesh2dx(&data_interp,index,x_data,y_data,nods_data,nels_data,data,data_rows,data_cols,x_interp,y_interp,nods_interp,default_values,num_default_values,contours,numcontours);
 
         /*Write data: */
@@ -88 +143 @@
         _printf_("\n");
         _printf_("   Usage:\n");
-        _printf_("      data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp,default_value);\n");
+        _printf_("         data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp);\n");
+        _printf_("      or data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp,default_value,contourname);\n");
         _printf_("\n");
         _printf_("      index: index of the mesh where data is defined\n");
@@ -94 +150 @@
         _printf_("      data: matrix holding the data to be interpolated onto the mesh. (one column per field)\n");
         _printf_("      x_interp,y_interp: coordinates of the points onto which we interpolate.\n");
-        _printf_("      default_value: default value if no data is found (holes).\n");
+        _printf_("      if default_value and contourname not specified: linear interpolation will happen on all x_interp,y_interp.\n");
+        _printf_("      if (default_value,contourname) specified: linear interpolation will happen on all x_interp,y_interp inside the contour, default value will be adopted on the rest of the mesh.\n");
+        _printf_("      note that default_value is either a scalar, or a vector of size  length(x_interp)\n");
         _printf_("      data_interp: vector of mesh interpolated data.\n");
         _printf_("\n");

issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.h ¶

@@ -25 +25 @@
 #define YINTERPHANDLE prhs[5]
 #define DEFAULTHANDLE prhs[6]
+#define FILENAME prhs[7]
 
 /* serial output macros: */
@@ -32 +33 @@
 #undef NLHS
 #define NLHS  1
-#undef NRHS
-#define NRHS  7
 
 #endif

issm/trunk/src/mex/Makefile.am ¶

@@ -35 +35 @@
                                 InterpFromMeshToMesh3d \
                                 InterpFromMeshToGrid \
+                                InterpFromMesh2d \
                                 MassFlux\
                                 Mergesolutionfromftog\
@@ -191 +192 @@
                                                                         InterpFromMeshToGrid/InterpFromMeshToGrid.h
 
+InterpFromMesh2d_SOURCES = InterpFromMesh2d/InterpFromMesh2d.cpp\
+                                                                        InterpFromMesh2d/InterpFromMesh2d.h
+
 AverageFilter_SOURCES = AverageFilter/AverageFilter.cpp\
                           AverageFilter/AverageFilter.h