source: issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp@ 8306

/*!\file InterpFromMesh2d.c
 * \brief: data interpolation from a list of (x,y,values) into mesh vertices
 
        InterpFromMesh2d.c

        usage:
        data_mesh=InterpFromMesh2d(index,x,y,data,x_mesh,y_mesh);
        
        where:

                input:
                x,y: coordinates of matrix data
                data - matrix holding the data to be interpolated onto the mesh.
                x_mesh,y_mesh: coordinates of the mesh vertices onto which we interpolate.
                
                output: 
                data_mesh:  vector of mesh interpolated data.

*/
        
#include "./InterpFromMesh2d.h"

void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {

        /*input: */
        double* index_data=NULL;
        int     index_data_rows;
        int     dummy;

        double* x_data=NULL;
        int     x_data_rows;
        
        double* y_data=NULL;
        int     y_data_rows;

        double* data=NULL; 
        int     data_rows;
        int     data_cols;

        double* x_prime=NULL;
        double* y_prime=NULL;
        
        int     x_prime_rows;
        int     y_prime_rows;


        double* default_values=NULL;
        int     num_default_values=0;

        //contours
        mxArray*  matlabstructure=NULL;
        Contour** contours=NULL;
        int       numcontours;
        Contour*  contouri=NULL;
        int       i;

        /*Intermediary*/
        int nods_data;
        int nels_data;
        int nods_prime;

        /* output: */
        Vec  data_prime=NULL;

        /*Boot module: */
        MODULEBOOT();

        /*checks on arguments on the matlab side: */
        if(nlhs!=NLHS){
                InterpFromMesh2dUsage();
                _error_("InterpFromMeshToMesh2dUsage usage error");
        }
        if((nrhs!=6) && (nrhs!=7) && (nrhs!=8)){
                InterpFromMesh2dUsage();
                _error_("InterpFromMeshToMesh2dUsage usage error");
        }

        /*Input datasets: */
        FetchData(&index_data,&index_data_rows,&dummy,INDEXHANDLE);
        FetchData(&x_data,&x_data_rows,NULL,XHANDLE);
        FetchData(&y_data,&y_data_rows,NULL,YHANDLE);
        FetchData(&data,&data_rows,&data_cols,DATAHANDLE);
        FetchData(&x_prime,&x_prime_rows,NULL,XPRIMEHANDLE);
        FetchData(&y_prime,&y_prime_rows,NULL,YPRIMEHANDLE);

        if(nrhs>=7){
                /*default values: */
                FetchData(&default_values,&num_default_values,DEFAULTHANDLE);
        }
        else{
                default_values=NULL;
                num_default_values=0;
        }

        if(nrhs>=8){
                
                /*Call expread on filename to build a contour array in the matlab workspace: */
                mexCallMATLAB( 1, &matlabstructure, 1, (mxArray**)&FILENAME, "expread");

                /*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 vertices %i\n",contouri->nods);
                  for (j=0;j<contouri->nods;j++){
                  printf("   %lf %lf\n",*(contouri->x+j),*(contouri->y+j));
                  }
                  }*/
                /*}}}*/
        }
        else{
                numcontours=0;
                contours=NULL;
        }


        /*some checks*/
        if (x_data_rows!=y_data_rows){
                _error_("vectors x and y should have the same length!");
        }
        if (x_prime_rows!=y_prime_rows){
                _error_("vectors x_prime and y_prime should have the same length!");
        }
        
        /*get number of elements and number of nodes in the data*/
        nels_data=index_data_rows;
        nods_data=x_data_rows;
        nods_prime=x_prime_rows;

        /* Run core computations: */
        InterpFromMesh2dx(&data_prime,index_data,x_data,y_data,nods_data,nels_data,data,data_rows,x_prime,y_prime,nods_prime,default_values,num_default_values,contours,numcontours);

        /*Write data: */
        WriteData(DATAPRIME,data_prime);

        /*end module: */
        MODULEEND();
}

void InterpFromMesh2dUsage(void)
{
        _printf_(true,"   usage:\n");
        _printf_(true,"         data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime);\n\n");
        _printf_(true,"      or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value);\n\n");
        _printf_(true,"      or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value,contourname);\n\n");
        _printf_(true,"   where:\n");
        _printf_(true,"      x,y: coordinates of the nodes where data is defined\n");
        _printf_(true,"      index: index of the mesh where data is defined\n");
        _printf_(true,"      data - vector holding the data to be interpolated onto the points.\n");
        _printf_(true,"      x_prime,y_prime: coordinates of the mesh vertices onto which we interpolate.\n");
        _printf_(true,"      default_value: a scalar or vector of size length(x_prime).\n");
        _printf_(true,"      contourname: 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_(true,"      data_prime:  vector of prime interpolated data.\n");
        _printf_(true,"\n");
}