Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 5031) +++ /issm/trunk/src/c/Makefile.am (revision 5032) @@ -479,4 +479,6 @@ ./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,4 +1017,6 @@ ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\ ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\ + ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\ + ./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\ ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\ ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\ Index: /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp =================================================================== --- /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp (revision 5032) +++ /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp (revision 5032) @@ -0,0 +1,138 @@ +/*!\file: InterpFromMesh2dx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#include "./InterpFromMesh2dx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../objects/objects.h" +#include "../modules.h" + +int InterpFromMesh2dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime, + double* default_values,int num_default_values,Contour** contours,int numcontours){ + + + /*Output*/ + Vec data_prime=NULL; + + /*Intermediary*/ + int i,j; + int interpolation_type; + bool debug; + double area; + double area_1,area_2,area_3; + double data_value; + double xmin,xmax; + double ymin,ymax; + + /*contours: */ + Vec vec_incontour=NULL; + double* incontour=NULL; + + /*some checks*/ + if (nels_data<1 || nods_data<3 || nods_prime==0){ + ISSMERROR("nothing to be done according to the mesh given in input"); + } + + /*Set debug to 1 if there are lots of elements*/ + debug=(bool)((double)nels_data*(double)nods_prime >= pow((double)10,(double)9)); + + /*figure out what kind of interpolation is needed*/ + if (data_length==nods_data){ + interpolation_type=1; + } + else if (data_length==nels_data){ + interpolation_type=2; + } + else{ + ISSMERROR("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!"); + } + + if((numcontours) && (interpolation_type==2)){ + ISSMERROR(" element interpolation_type with contours not supported yet!"); + } + + /*Get prime mesh extrema coordinates*/ + xmin=x_prime[0]; xmax=x_prime[0];ymin=y_prime[0]; ymax=y_prime[0]; + for (i=1;ixmax) xmax=x_prime[i]; + if (y_prime[i]ymax) ymax=y_prime[i]; + } + + /*Initialize output*/ + data_prime=NewVec(nods_prime); + if(num_default_values){ + if(num_default_values==1)for (i=0;ixmax) && (x_data[(int)index_data[3*i+1]-1]>xmax) && (x_data[(int)index_data[3*i+2]-1]>xmax)) continue; + if ( (y_data[(int)index_data[3*i+0]-1]ymax) && (y_data[(int)index_data[3*i+1]-1]>ymax) && (y_data[(int)index_data[3*i+2]-1]>ymax)) continue; + + /*get area of the current element (Jacobian = 2 * area)*/ + //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1; + area=x_data[(int)index_data[3*i+1]-1]*y_data[(int)index_data[3*i+2]-1]-y_data[(int)index_data[3*i+1]-1]*x_data[(int)index_data[3*i+2]-1] + + x_data[(int)index_data[3*i+0]-1]*y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+0]-1]*x_data[(int)index_data[3*i+1]-1] + + x_data[(int)index_data[3*i+2]-1]*y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1]*x_data[(int)index_data[3*i+0]-1]; + + /*loop over the prime nodes*/ + for (j=0;j=0 && area_2>=0 && area_3>=0){ + + /*Yes ! compute the value on the point*/ + if (interpolation_type==1){ + /*nodal interpolation*/ + data_value=area_1*data[(int)index_data[3*i+0]-1]+area_2*data[(int)index_data[3*i+1]-1]+area_3*data[(int)index_data[3*i+2]-1]; + } + else{ + /*element interpolation*/ + data_value=data[i]; + } + if (isnan(data_value)){ + if(num_default_values==1) data_value=default_values[0]; + else data_value=default_values[j]; + } + + /*insert value and go to the next point*/ + VecSetValue(data_prime,j,data_value,INSERT_VALUES); + } + } + } + } + if (debug) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0); + + /*Assign output pointers:*/ + *pdata_prime=data_prime; +} Index: /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h =================================================================== --- /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h (revision 5032) +++ /issm/trunk/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dx.h (revision 5032) @@ -0,0 +1,15 @@ +/*!\file InterpFromMesh2dx.h + * \brief: header file for Data interpolation routines. + */ + +#ifndef _INTERPFROMMESH2DX_H +#define _INTERPFROMMESH2DX_H + +#include "../../objects/objects.h" +#include "../../toolkits/toolkits.h" + +int InterpFromMesh2dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, int nods_prime, + double* default_values,int num_default_values,Contour** contours,int numcontours); + +#endif /* _INTERPFROMMESH2DX_H */ + Index: /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp =================================================================== --- /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (revision 5031) +++ /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (revision 5032) @@ -8,4 +8,5 @@ #include "../../toolkits/toolkits.h" #include "../../objects/objects.h" +#include "../modules.h" using namespace bamg; @@ -13,6 +14,6 @@ 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,4 +31,10 @@ int verbose=0; + /*default values: */ + Vec vec_incontour=NULL; + double* incontour=NULL; + bool skip_bamg=false; + + /*Checks*/ if (data_cols<=0){ @@ -36,4 +43,13 @@ 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,59 +66,74 @@ if (verbose) printf("Loop over the nodes\n"); for(i=0;i0){ - //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;j0){ + //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=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=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]; } - } Index: /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h =================================================================== --- /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h (revision 5031) +++ /issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h (revision 5032) @@ -10,5 +10,5 @@ /* 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 Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 5031) +++ /issm/trunk/src/c/modules/modules.h (revision 5032) @@ -30,4 +30,5 @@ #include "./InputDuplicatex/InputDuplicatex.h" #include "./InputScalex/InputScalex.h" +#include "./InterpFromMesh2dx/InterpFromMesh2dx.h" #include "./InterpFromGridToMeshx/InterpFromGridToMeshx.h" #include "./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h" Index: /issm/trunk/src/mex/InterpFromMesh2d/1 =================================================================== --- /issm/trunk/src/mex/InterpFromMesh2d/1 (revision 5032) +++ /issm/trunk/src/mex/InterpFromMesh2d/1 (revision 5032) @@ -0,0 +1,169 @@ +/*!\file InterpFromMesh2d.c + * \brief: data interpolation from a list of (x,y,values) into mesh grids + + 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 grids 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(); + ISSMERROR("InterpFromMeshToMesh2dUsage usage error"); + } + if((nrhs!=6) && (nrhs!=7) && (nrhs!=8)){ + InterpFromMesh2dUsage(); + ISSMERROR("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;inods=(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;inods); + for (j=0;jnods;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){ + ISSMERROR("vectors x and y should have the same length!"); + } + if (x_prime_rows!=y_prime_rows){ + ISSMERROR("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_(" usage:\n"); + _printf_(" data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime);\n\n"); + _printf_(" or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value);\n\n"); + _printf_(" or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value,contourname);\n\n"); + _printf_(" where:\n"); + _printf_(" x,y: coordinates of the nodes where data is defined\n"); + _printf_(" index: index of the mesh where data is defined\n"); + _printf_(" data - vector holding the data to be interpolated onto the points.\n"); + _printf_(" x_prime,y_prime: coordinates of the mesh grids onto which we interpolate.\n"); + _printf_(" default_value: a scalar or vector of size length(x_prime).\n"); + _printf_(" 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_(" data_prime: vector of prime interpolated data.\n"); + _printf_("\n"); +} Index: /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp =================================================================== --- /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp (revision 5032) +++ /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp (revision 5032) @@ -0,0 +1,169 @@ +/*!\file InterpFromMesh2d.c + * \brief: data interpolation from a list of (x,y,values) into mesh grids + + 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 grids 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(); + ISSMERROR("InterpFromMeshToMesh2dUsage usage error"); + } + if((nrhs!=6) && (nrhs!=7) && (nrhs!=8)){ + InterpFromMesh2dUsage(); + ISSMERROR("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;inods=(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;inods); + for (j=0;jnods;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){ + ISSMERROR("vectors x and y should have the same length!"); + } + if (x_prime_rows!=y_prime_rows){ + ISSMERROR("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_(" usage:\n"); + _printf_(" data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime);\n\n"); + _printf_(" or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value);\n\n"); + _printf_(" or data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime,default_value,contourname);\n\n"); + _printf_(" where:\n"); + _printf_(" x,y: coordinates of the nodes where data is defined\n"); + _printf_(" index: index of the mesh where data is defined\n"); + _printf_(" data - vector holding the data to be interpolated onto the points.\n"); + _printf_(" x_prime,y_prime: coordinates of the mesh grids onto which we interpolate.\n"); + _printf_(" default_value: a scalar or vector of size length(x_prime).\n"); + _printf_(" 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_(" data_prime: vector of prime interpolated data.\n"); + _printf_("\n"); +} Index: /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.h =================================================================== --- /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.h (revision 5032) +++ /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.h (revision 5032) @@ -0,0 +1,38 @@ +/*!\file InterpFromMesh2d.h + * \brief: prototype for Data Interpolation mex module. + */ + +#ifndef _INTERPFROMMESH2D_H +#define _INTERPFROMMESH2D_H + +/* local prototypes: */ +void InterpFromMesh2dUsage(void); + +#include "../../c/modules/modules.h" +#include "../../c/Container/Container.h" +#include "../../c/shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "InterpFromMesh2d" + +#undef CLEANUP +#define CLEANUP InterpFromMesh2dLocalCleanup(); + +/* serial input macros: */ +#define INDEXHANDLE prhs[0] +#define XHANDLE prhs[1] +#define YHANDLE prhs[2] +#define DATAHANDLE prhs[3] +#define XPRIMEHANDLE prhs[4] +#define YPRIMEHANDLE prhs[5] +#define DEFAULTHANDLE prhs[6] +#define FILENAME prhs[7] + +/* serial output macros: */ +#define DATAPRIME (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 + +#endif /* _INTERPFROMMESH2D_H */ Index: /issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp =================================================================== --- /issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp (revision 5031) +++ /issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp (revision 5032) @@ -5,4 +5,6 @@ void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + int i; /*input: */ @@ -26,5 +28,12 @@ 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,5 +50,12 @@ /*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,5 +67,44 @@ 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;inods=(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;inods); + for (j=0;jnods;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,5 +126,5 @@ /* 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,5 +143,6 @@ _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,5 +150,7 @@ _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"); Index: /issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.h =================================================================== --- /issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.h (revision 5031) +++ /issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.h (revision 5032) @@ -25,4 +25,5 @@ #define YINTERPHANDLE prhs[5] #define DEFAULTHANDLE prhs[6] +#define FILENAME prhs[7] /* serial output macros: */ @@ -32,6 +33,4 @@ #undef NLHS #define NLHS 1 -#undef NRHS -#define NRHS 7 #endif Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 5031) +++ /issm/trunk/src/mex/Makefile.am (revision 5032) @@ -35,4 +35,5 @@ InterpFromMeshToMesh3d \ InterpFromMeshToGrid \ + InterpFromMesh2d \ MassFlux\ Mergesolutionfromftog\ @@ -191,4 +192,7 @@ InterpFromMeshToGrid/InterpFromMeshToGrid.h +InterpFromMesh2d_SOURCES = InterpFromMesh2d/InterpFromMesh2d.cpp\ + InterpFromMesh2d/InterpFromMesh2d.h + AverageFilter_SOURCES = AverageFilter/AverageFilter.cpp\ AverageFilter/AverageFilter.h