Index: /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.cpp =================================================================== --- /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.cpp (revision 1172) +++ /issm/trunk/src/c/InterpFromGridx/InterpFromGridx.cpp (revision 1172) @@ -0,0 +1,183 @@ +/*!\file: InterpFromGridx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#include "./InterpFromGridx.h" +#include "../shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "InterpFromGridx" + +int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid); + +int InterpFromGridx( Vec* pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods) { + + + /*output: */ + Vec data_mesh=NULL; + + /*Intermediary*/ + double* x=NULL; + double* y=NULL; + int i,m,n; + int ind; + double dis; + double x_grid,y_grid; + double xi,eta; + double G1,G2,G3,G4,data_value; + double area_1,area_2,area_3; + double area; + + /*Some checks on arguments: */ + if ((M<=2) || (N<=2) || (nods<=0)){ + throw ErrorException(__FUNCT__,"nothing to be done according to the dimensions of input matrices and vectors."); + } + + /*Allocate output vector: */ + data_mesh=NewVec(nods); + + /*Find out what kind of coordinates (x_in,y_in) have been given is input*/ + if(N==(x_rows-1) && M==(y_rows-1)){ + + /*The coordinates given in input describe the contour of each pixel. Take the center of each pixel*/ + x=(double*)xmalloc(N*sizeof(double)); + y=(double*)xmalloc(M*sizeof(double)); + for (i=0;ixmax) xmax=x_prime[i]; + if (y_prime[i]ymax) ymax=y_prime[i]; + } + + /*Initialize output*/ + data_prime=NewVec(nods_prime); + + printf("\n interpolation progress: %5.2lf %%",0.0); + 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;j0 || area_1==0) && (area_2>0 || area_2==0) && (area_3>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]; + } + + /*insert value and go to the next point*/ + VecSetValue(data_prime,j,data_value,INSERT_VALUES); + } + } + } + printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0); + + /*Write data: */ + WriteData(DATAPRIME,data_prime,0,0,"Vector",NULL); + + /*end module: */ + MODULEEND(); +} + +void InterpFromMesh2dUsage(void) +{ + _printf_(" usage:\n"); + _printf_(" data_prime=InterpFromMesh2d(index,x,y,data,x_prime,y_prime);\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_(" 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 1172) +++ /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.h (revision 1172) @@ -0,0 +1,37 @@ +/*!\file InterpFromMesh2d.h + * \brief: prototype for Data Interpolation mex module. + */ + +#ifndef _INTERPFROMMESH2D_H +#define _INTERPFROMMESH2D_H + +/* local prototypes: */ +void InterpFromMesh2dUsage(void); + +#include "../../c/issm.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] + +/* serial output macros: */ +#define DATAPRIME (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 +#undef NRHS +#define NRHS 6 + +#endif /* _DATAINTERP2_H */ Index: /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp =================================================================== --- /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp (revision 1172) +++ /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp (revision 1172) @@ -0,0 +1,195 @@ +/*!\file InterpFromMesh3d.c + * \brief: data interpolation from a list of (x,y,values) into mesh grids + + InterpFromMesh3d.c + + usage: + data_mesh=InterpFromMesh3d(index,x,y,z,data,x_mesh,y_mesh,z_mesh); + + where: + + input: + x,y,z: coordinates of matrix data + data - matrix holding the data to be interpolated onto the mesh. + x_mesh,y_mesh,z_mesh: coordinates of the mesh grids onto which we interpolate. + + output: + data_mesh: vector of mesh interpolated data. + +*/ + +#include "./InterpFromMesh3d.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + + /*input: */ + double* index_data=NULL; + int index_data_rows; + + double* x_data=NULL; + double* y_data=NULL; + double* z_data=NULL; + + int x_data_rows; + int y_data_rows; + int z_data_rows; + + double* data=NULL; + int data_rows; + int data_cols; + + double* x_prime=NULL; + double* y_prime=NULL; + double* z_prime=NULL; + + int x_prime_rows; + int y_prime_rows; + int z_prime_rows; + + /*Intermediary*/ + int i,j; + int nods_data; + int nels_data; + int nods_prime; + int interpolation_type; + double area; + double area_1,area_2,area_3; + double zeta,bed,surface; + double data_value; + double xmin,xmax; + double ymin,ymax; + + /* output: */ + Vec data_prime=NULL; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InterpFromMesh3dUsage); + + /*Input datasets: */ + FetchData((void**)&index_data,&index_data_rows,NULL,INDEXHANDLE,"Matrix","Mat"); + FetchData((void**)&x_data,&x_data_rows,NULL,XHANDLE,"Matrix","Mat"); + FetchData((void**)&y_data,&y_data_rows,NULL,YHANDLE,"Matrix","Mat"); + FetchData((void**)&z_data,&z_data_rows,NULL,ZHANDLE,"Matrix","Mat"); + FetchData((void**)&data,&data_rows,&data_cols,DATAHANDLE,"Matrix","Mat"); + FetchData((void**)&x_prime,&x_prime_rows,NULL,XPRIMEHANDLE,"Matrix","Mat"); + FetchData((void**)&y_prime,&y_prime_rows,NULL,YPRIMEHANDLE,"Matrix","Mat"); + FetchData((void**)&z_prime,&z_prime_rows,NULL,ZPRIMEHANDLE,"Matrix","Mat"); + + /*some checks*/ + if (index_data_rows<1 || x_data_rows<6 || y_data_rows<6 || z_data_rows<6){ + throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input"); + } + if (x_data_rows!=y_data_rows || x_data_rows!=z_data_rows){ + throw ErrorException(__FUNCT__,"vectors x, y and z should have the same length!"); + } + if (x_prime_rows!=y_prime_rows || x_prime_rows!=z_prime_rows){ + throw ErrorException(__FUNCT__,"vectors x_prime, y_prime and z_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; + + /*figure out what kind of interpolation is needed*/ + if (data_rows==nods_data){ + interpolation_type=1; + } + else if (data_rows==nels_data){ + interpolation_type=2; + } + else{ + throw ErrorException(__FUNCT__,"length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!"); + } + + /*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); + + printf("\n interpolation progress: %5.2lf %%",0.0); + for (i=0;ixmax) && (x_data[(int)index_data[6*i+1]-1]>xmax) && (x_data[(int)index_data[6*i+2]-1]>xmax)) continue; + if ( (y_data[(int)index_data[6*i+0]-1]ymax) && (y_data[(int)index_data[6*i+1]-1]>ymax) && (y_data[(int)index_data[6*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[6*i+1]-1]*y_data[(int)index_data[6*i+2]-1]-y_data[(int)index_data[6*i+1]-1]*x_data[(int)index_data[6*i+2]-1]+ x_data[(int)index_data[6*i+0]-1]*y_data[(int)index_data[6*i+1]-1]-y_data[(int)index_data[6*i+0]-1]*x_data[(int)index_data[6*i+1]-1]+ x_data[(int)index_data[6*i+2]-1]*y_data[(int)index_data[6*i+0]-1]-y_data[(int)index_data[6*i+2]-1]*x_data[(int)index_data[6*i+0]-1]; + + /*loop over the prime nodes*/ + for (j=0;j0 || area_1==0) && (area_2>0 || area_2==0) && (area_3>0 || area_3==0)){ + + /*compute bottom and top height of the element at this 2d position*/ + bed =area_1*z_data[(int)index_data[6*i+0]-1]+area_2*z_data[(int)index_data[6*i+1]-1]+area_3*z_data[(int)index_data[6*i+2]-1]; + surface=area_1*z_data[(int)index_data[6*i+3]-1]+area_2*z_data[(int)index_data[6*i+4]-1]+area_3*z_data[(int)index_data[6*i+5]-1]; + + /*Compute zeta*/ + zeta=2*(z_prime[j]-bed)/(surface-bed)-1; + + if (zeta >=-1 && zeta<=1){ + if (interpolation_type==1){ + /*nodal interpolation*/ + data_value=(1-zeta)/2*(area_1*data[(int)index_data[6*i+0]-1]+area_2*data[(int)index_data[6*i+1]-1]+area_3*data[(int)index_data[6*i+2]-1]) + (1+zeta)/2*(area_1*data[(int)index_data[6*i+3]-1]+area_2*data[(int)index_data[6*i+4]-1]+area_3*data[(int)index_data[6*i+5]-1]); + } + else{ + /*element interpolation*/ + data_value=data[i]; + } + + /*insert value and go to the next point*/ + VecSetValue(data_prime,j,data_value,INSERT_VALUES); + } + } + } + } + printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0); + + /*Write data: */ + WriteData(DATAPRIME,data_prime,0,0,"Vector",NULL); + + /*end module: */ + MODULEEND(); +} + +void InterpFromMesh3dUsage(void) +{ + _printf_(" usage:\n"); + _printf_(" data_prime=InterpFromMesh3d(index,x,y,z,data,x_prime,y_prime,z_prime);\n\n"); + _printf_(" where:\n"); + _printf_(" x,y,z: 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,z_prime: coordinates of the mesh grids onto which we interpolate.\n"); + _printf_(" data_prime: vector of prime interpolated data.\n"); + _printf_("\n"); +} Index: /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.h =================================================================== --- /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.h (revision 1172) +++ /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.h (revision 1172) @@ -0,0 +1,39 @@ +/*!\file InterpFromMesh3d.h + * \brief: prototype for Data Interpolation mex module. + */ + +#ifndef _INTERPFROMMESH3D_H +#define _INTERPFROMMESH3D_H + +/* local prototypes: */ +void InterpFromMesh3dUsage(void); + +#include "../../c/issm.h" + +#undef __FUNCT__ +#define __FUNCT__ "InterpFromMesh3d" + +#undef CLEANUP +#define CLEANUP InterpFromMesh3dLocalCleanup(); + + +/* serial input macros: */ +#define INDEXHANDLE prhs[0] +#define XHANDLE prhs[1] +#define YHANDLE prhs[2] +#define ZHANDLE prhs[3] +#define DATAHANDLE prhs[4] +#define XPRIMEHANDLE prhs[5] +#define YPRIMEHANDLE prhs[6] +#define ZPRIMEHANDLE prhs[7] + +/* serial output macros: */ +#define DATAPRIME (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 +#undef NRHS +#define NRHS 8 + +#endif /* _INTERPFROMMESH3D_H */ Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 1171) +++ /issm/trunk/src/mex/Makefile.am (revision 1172) @@ -12,5 +12,4 @@ ContourToNodes \ ControlConstrain \ - DataInterp \ Dof\ Du\ @@ -20,4 +19,7 @@ GriddataMeshToGrid\ HoleFiller \ + InterpFromGrid \ + InterpFromMesh2d \ + InterpFromMesh3d \ Mergesolutionfromftog\ MeshPartition\ @@ -96,7 +98,4 @@ ControlConstrain/ControlConstrain.h -DataInterp_SOURCES = DataInterp/DataInterp.cpp\ - DataInterp/DataInterp.h - Dof_SOURCES = Dof/Dof.cpp\ Dof/Dof.h @@ -123,4 +122,13 @@ HoleFiller/HoleFiller.h +InterpFromGrid_SOURCES = InterpFromGrid/InterpFromGrid.cpp\ + InterpFromGrid/InterpFromGrid.h + +InterpFromMesh2d_SOURCES = InterpFromMesh2d/InterpFromMesh2d.cpp\ + InterpFromMesh2d/InterpFromMesh2d.h + +InterpFromMesh3d_SOURCES = InterpFromMesh3d/InterpFromMesh3d.cpp\ + InterpFromMesh3d/InterpFromMesh3d.h + Mergesolutionfromftog_SOURCES = Mergesolutionfromftog/Mergesolutionfromftog.cpp\ Mergesolutionfromftog/Mergesolutionfromftog.h