Index: sm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dToMeshx.cpp =================================================================== --- /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dToMeshx.cpp (revision 2291) +++ (revision ) @@ -1,112 +1,0 @@ -/*!\file: InterpFromMeshToMesh2dx.cpp - * \brief "c" core code for interpolating values from a structured grid. - */ - -#include "./InterpFromMeshToMesh2dx.h" -#include "../shared/shared.h" - -#undef __FUNCT__ -#define __FUNCT__ "InterpFromMeshToMesh2dx" - -int InterpFromMeshToMesh2dx( 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; - bool debug; - 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"); - } - - /*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{ - 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); - } - } - } - if (debug) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0); - - /*Assign output pointers:*/ - *pdata_prime=data_prime; -} Index: sm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dToMeshx.h =================================================================== --- /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMesh2dToMeshx.h (revision 2291) +++ (revision ) @@ -1,13 +1,0 @@ -/*!\file InterpFromMeshToMesh2dx.h - * \brief: header file for Data interpolation routines. - */ - -#ifndef _INTERPFROMMESHTOMESH2DX_H -#define _INTERPFROMMESHTOMESH2DX_H - -#include "../toolkits/toolkits.h" - -int InterpFromMeshToMesh2dx( 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 /* _INTERPFROMMESHTOMESH2DX_H */ - Index: /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp =================================================================== --- /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (revision 2292) +++ /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (revision 2292) @@ -0,0 +1,112 @@ +/*!\file: InterpFromMeshToMesh2dx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#include "./InterpFromMeshToMesh2dx.h" +#include "../shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "InterpFromMeshToMesh2dx" + +int InterpFromMeshToMesh2dx( 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; + bool debug; + 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"); + } + + /*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{ + 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); + } + } + } + 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/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h =================================================================== --- /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h (revision 2292) +++ /issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h (revision 2292) @@ -0,0 +1,13 @@ +/*!\file InterpFromMeshToMesh2dx.h + * \brief: header file for Data interpolation routines. + */ + +#ifndef _INTERPFROMMESHTOMESH2DX_H +#define _INTERPFROMMESHTOMESH2DX_H + +#include "../toolkits/toolkits.h" + +int InterpFromMeshToMesh2dx( 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 /* _INTERPFROMMESHTOMESH2DX_H */ + Index: sm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dToMeshx.cpp =================================================================== --- /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dToMeshx.cpp (revision 2291) +++ (revision ) @@ -1,122 +1,0 @@ -/*!\file: InterpFromMeshToMesh3dx.cpp - * \brief "c" core code for interpolating values from a structured grid. - */ - -#include "./InterpFromMeshToMesh3dx.h" -#include "../shared/shared.h" - -#undef __FUNCT__ -#define __FUNCT__ "InterpFromMeshToMesh3dx" - -int InterpFromMeshToMesh3dx( 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; - bool debug; - 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"); - } - - /*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{ - 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); - } - } - } - } - if (debug) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0); - - /*Assign output pointers:*/ - *pdata_prime=data_prime; -} Index: sm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dToMeshx.h =================================================================== --- /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMesh3dToMeshx.h (revision 2291) +++ (revision ) @@ -1,13 +1,0 @@ -/*!\file InterpFromMeshToMesh3dx.h - * \brief: header file for Data interpolation routines. - */ - -#ifndef _INTERPFROMMESHTOMESH3DX_H -#define _INTERPFROMMESHTOMESH3DX_H - -#include "../toolkits/toolkits.h" - -int InterpFromMeshToMesh3dx( 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 /* _INTERPFROMMESHTOMESH3DX_H */ - Index: /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp =================================================================== --- /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp (revision 2292) +++ /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp (revision 2292) @@ -0,0 +1,122 @@ +/*!\file: InterpFromMeshToMesh3dx.cpp + * \brief "c" core code for interpolating values from a structured grid. + */ + +#include "./InterpFromMeshToMesh3dx.h" +#include "../shared/shared.h" + +#undef __FUNCT__ +#define __FUNCT__ "InterpFromMeshToMesh3dx" + +int InterpFromMeshToMesh3dx( 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; + bool debug; + 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"); + } + + /*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{ + 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); + } + } + } + } + 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/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h =================================================================== --- /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h (revision 2292) +++ /issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.h (revision 2292) @@ -0,0 +1,13 @@ +/*!\file InterpFromMeshToMesh3dx.h + * \brief: header file for Data interpolation routines. + */ + +#ifndef _INTERPFROMMESHTOMESH3DX_H +#define _INTERPFROMMESHTOMESH3DX_H + +#include "../toolkits/toolkits.h" + +int InterpFromMeshToMesh3dx( 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 /* _INTERPFROMMESHTOMESH3DX_H */ +