Index: /issm/trunk/src/c/InterpFromMesh2dx/InterpFromMesh2dx.cpp =================================================================== --- /issm/trunk/src/c/InterpFromMesh2dx/InterpFromMesh2dx.cpp (revision 1189) +++ /issm/trunk/src/c/InterpFromMesh2dx/InterpFromMesh2dx.cpp (revision 1189) @@ -0,0 +1,110 @@ +/*!\file: InterpFromMesh2dx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#include "./InterpFromMesh2dx.h" +#include "../shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "InterpFromMesh2dx" + +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_value) { + + /*Output*/ + Vec data_prime=NULL; + + /*Intermediary*/ + int i,j; + int interpolation_type; + double area; + double area_1,area_2,area_3; + double data_value; + double xmin,xmax; + double ymin,ymax; + + /*some checks*/ + if (nels_data<1 || nods_data<3 || nods_prime==0){ + throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input"); + } + + /*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{ + 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); + 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) data_value=default_value; + + /*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); + + /*Assign output pointers:*/ + *pdata_prime=data_prime; +} Index: /issm/trunk/src/c/InterpFromMesh2dx/InterpFromMesh2dx.h =================================================================== --- /issm/trunk/src/c/InterpFromMesh2dx/InterpFromMesh2dx.h (revision 1189) +++ /issm/trunk/src/c/InterpFromMesh2dx/InterpFromMesh2dx.h (revision 1189) @@ -0,0 +1,13 @@ +/*!\file InterpFromMesh2dx.h + * \brief: header file for Data interpolation routines. + */ + +#ifndef _INTERPFROMMESH2DX_H +#define _INTERPFROMMESH2DX_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_value); + +#endif /* _INTERPFROMMESH2DX_H */ + Index: /issm/trunk/src/c/InterpFromMesh3dx/InterpFromMesh3dx.cpp =================================================================== --- /issm/trunk/src/c/InterpFromMesh3dx/InterpFromMesh3dx.cpp (revision 1189) +++ /issm/trunk/src/c/InterpFromMesh3dx/InterpFromMesh3dx.cpp (revision 1189) @@ -0,0 +1,120 @@ +/*!\file: InterpFromMesh3dx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#include "./InterpFromMesh3dx.h" +#include "../shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "InterpFromMesh3dx" + +int InterpFromMesh3dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value) { + + /*Output*/ + Vec data_prime=NULL; + + /*Intermediary*/ + int i,j; + 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; + + /*some checks*/ + if (nels_data<1 || nods_data<6 || nods_prime==0){ + throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input"); + } + + /*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{ + 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); + 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;j=0 && area_2>=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]; + } + if isnan(data_value) data_value=default_value; + + /*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); + + /*Assign output pointers:*/ + *pdata_prime=data_prime; +} Index: /issm/trunk/src/c/InterpFromMesh3dx/InterpFromMesh3dx.h =================================================================== --- /issm/trunk/src/c/InterpFromMesh3dx/InterpFromMesh3dx.h (revision 1189) +++ /issm/trunk/src/c/InterpFromMesh3dx/InterpFromMesh3dx.h (revision 1189) @@ -0,0 +1,13 @@ +/*!\file InterpFromMesh3dx.h + * \brief: header file for Data interpolation routines. + */ + +#ifndef _INTERPFROMMESH3DX_H +#define _INTERPFROMMESH3DX_H + +#include "../toolkits/toolkits.h" + +int InterpFromMesh3dx( Vec* pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value); + +#endif /* _INTERPFROMMESH3DX_H */ + Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 1188) +++ /issm/trunk/src/c/Makefile.am (revision 1189) @@ -241,4 +241,8 @@ ./InterpFromGridx/InterpFromGridx.cpp\ ./InterpFromGridx/InterpFromGridx.h\ + ./InterpFromMesh2dx/InterpFromMesh2dx.cpp\ + ./InterpFromMesh2dx/InterpFromMesh2dx.h\ + ./InterpFromMesh3dx/InterpFromMesh3dx.cpp\ + ./InterpFromMesh3dx/InterpFromMesh3dx.h\ ./HoleFillerx/HoleFillerx.cpp\ ./HoleFillerx/HoleFillerx.h\ @@ -521,4 +525,8 @@ ./InterpFromGridx/InterpFromGridx.cpp\ ./InterpFromGridx/InterpFromGridx.h\ + ./InterpFromMesh2dx/InterpFromMesh2dx.cpp\ + ./InterpFromMesh2dx/InterpFromMesh2dx.h\ + ./InterpFromMesh3dx/InterpFromMesh3dx.cpp\ + ./InterpFromMesh3dx/InterpFromMesh3dx.h\ ./HoleFillerx/HoleFillerx.cpp\ ./HoleFillerx/HoleFillerx.h\ Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 1188) +++ /issm/trunk/src/c/issm.h (revision 1189) @@ -26,4 +26,6 @@ #include "./ContourToNodesx/ContourToNodesx.h" #include "./InterpFromGridx/InterpFromGridx.h" +#include "./InterpFromMesh2dx/InterpFromMesh2dx.h" +#include "./InterpFromMesh3dx/InterpFromMesh3dx.h" #include "./HoleFillerx/HoleFillerx.h" #include "./MeshPartitionx/MeshPartitionx.h" Index: /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp =================================================================== --- /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp (revision 1188) +++ /issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp (revision 1189) @@ -46,14 +46,7 @@ /*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 data_value; - double xmin,xmax; - double ymin,ymax; /* output: */ @@ -75,11 +68,5 @@ FetchData((void**)&default_value,NULL,NULL,DEFAULTHANDLE,"Scalar",NULL); - /* Run core computations: */ - //InterpFromMesh2dx( &data_mesh, index , index_rows, x, y, x_rows, data, data_rows, x_mesh, y_mesh, x_mesh_rows); - /*some checks*/ - if (index_data_rows<1 || x_data_rows<3 || y_data_rows<3){ - throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input"); - } if (x_data_rows!=y_data_rows){ throw ErrorException(__FUNCT__,"vectors x and y should have the same length!"); @@ -94,77 +81,6 @@ 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); - 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) data_value=default_value; - - /*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); + /* 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_value); /*Write data: */ Index: /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp =================================================================== --- /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp (revision 1188) +++ /issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp (revision 1189) @@ -50,15 +50,7 @@ /*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: */ @@ -83,7 +75,4 @@ /*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!"); @@ -92,5 +81,4 @@ 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; @@ -98,86 +86,6 @@ 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); - 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;j=0 && area_2>=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]; - } - if isnan(data_value) data_value=default_value; - - /*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); + /* Run core computations: */ + InterpFromMesh3dx(&data_prime,index_data,x_data,y_data,z_data,nods_data,nels_data,data,data_rows,x_prime,y_prime,z_prime,nods_prime,default_value); /*Write data: */