Index: sm/trunk/src/c/modules/TriaSearchx/InterpFromMeshToMesh2dx.cpp =================================================================== --- /issm/trunk/src/c/modules/TriaSearchx/InterpFromMeshToMesh2dx.cpp (revision 5350) +++ (revision ) @@ -1,146 +1,0 @@ -/*!\file InterpFromMeshToMesh2dx - */ - -#include "./InterpFromMeshToMesh2dx.h" - -#include "../../shared/shared.h" -#include "../../include/include.h" -#include "../../toolkits/toolkits.h" -#include "../../objects/objects.h" -#include "../modules.h" - -using namespace bamg; -using namespace std; - -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_values,int num_default_values, Contour** contours, int numcontours){ - - /*Output*/ - double* data_interp=NULL; - - /*Intermediary*/ - R2 r; - I2 I; - int i,j,k; - int it; - int i0,i1,i2; - double areacoord[3]; - double aa,bb; - double data_value; - Icoor2 dete[3]; - int verbose=0; - - /*default values: */ - Vec vec_incontour=NULL; - double* incontour=NULL; - bool skip_bamg=false; - - - /*Checks*/ - if (data_cols<=0){ - ISSMERROR("data provided has a negative number of columns"); - } - 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); - } - - /*Initialize output*/ - if (verbose) printf("Initializing output vector\n"); - data_interp=(double*)xmalloc(nods_interp*data_cols*sizeof(double)); - - // read background mesh - if (verbose) printf("Reading mesh\n"); - Mesh Th(index_data,x_data,y_data,nods_data,nels_data); - Th.CreateSingleVertexToTriangleConnectivity(); - - //Loop over output nodes - 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.GetId(tb[0]); - i1=Th.GetId(tb[1]); - i2=Th.GetId(tb[2]); - //triangle number - it=Th.GetId(tb); - } - //external point - else { - //Get closest adjacent triangle (inside the mesh) - AdjacentTriangle 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.GetId(tc[0]); - i1=Th.GetId(tc[1]); - i2=Th.GetId(tc[2]); - //triangle number - it=Th.GetId(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{ - if(num_default_values==1) data_interp[i]=default_values[0]; - else data_interp[i]=default_values[i]; - } - } - - /*Assign output pointers:*/ - if (verbose) printf("Assigning output\n"); - *pdata_interp=data_interp; - - /*No error return*/ - return 1; -} Index: sm/trunk/src/c/modules/TriaSearchx/InterpFromMeshToMesh2dx.h =================================================================== --- /issm/trunk/src/c/modules/TriaSearchx/InterpFromMeshToMesh2dx.h (revision 5350) +++ (revision ) @@ -1,14 +1,0 @@ -/*!\file: InterpFromMeshToMesh2dx.h - * \brief header file for Bamg module - */ - -#ifndef _INTERPFROMMESHTOMESH2DX_H -#define _INTERPFROMMESHTOMESH2DX_H - -#include "../../objects/objects.h" - -/* 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_values,int num_default_values,Contour** contours, int numcontours ); - -#endif Index: /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.cpp =================================================================== --- /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.cpp (revision 5351) +++ /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.cpp (revision 5351) @@ -0,0 +1,146 @@ +/*!\file InterpFromMeshToMesh2dx + */ + +#include "./InterpFromMeshToMesh2dx.h" + +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../objects/objects.h" +#include "../modules.h" + +using namespace bamg; +using namespace std; + +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_values,int num_default_values, Contour** contours, int numcontours){ + + /*Output*/ + double* data_interp=NULL; + + /*Intermediary*/ + R2 r; + I2 I; + int i,j,k; + int it; + int i0,i1,i2; + double areacoord[3]; + double aa,bb; + double data_value; + Icoor2 dete[3]; + int verbose=0; + + /*default values: */ + Vec vec_incontour=NULL; + double* incontour=NULL; + bool skip_bamg=false; + + + /*Checks*/ + if (data_cols<=0){ + ISSMERROR("data provided has a negative number of columns"); + } + 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); + } + + /*Initialize output*/ + if (verbose) printf("Initializing output vector\n"); + data_interp=(double*)xmalloc(nods_interp*data_cols*sizeof(double)); + + // read background mesh + if (verbose) printf("Reading mesh\n"); + Mesh Th(index_data,x_data,y_data,nods_data,nels_data); + Th.CreateSingleVertexToTriangleConnectivity(); + + //Loop over output nodes + 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.GetId(tb[0]); + i1=Th.GetId(tb[1]); + i2=Th.GetId(tb[2]); + //triangle number + it=Th.GetId(tb); + } + //external point + else { + //Get closest adjacent triangle (inside the mesh) + AdjacentTriangle 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.GetId(tc[0]); + i1=Th.GetId(tc[1]); + i2=Th.GetId(tc[2]); + //triangle number + it=Th.GetId(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{ + if(num_default_values==1) data_interp[i]=default_values[0]; + else data_interp[i]=default_values[i]; + } + } + + /*Assign output pointers:*/ + if (verbose) printf("Assigning output\n"); + *pdata_interp=data_interp; + + /*No error return*/ + return 1; +} Index: /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.h =================================================================== --- /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.h (revision 5351) +++ /issm/trunk/src/c/modules/TriaSearchx/TriaSearchx.h (revision 5351) @@ -0,0 +1,14 @@ +/*!\file: InterpFromMeshToMesh2dx.h + * \brief header file for Bamg module + */ + +#ifndef _INTERPFROMMESHTOMESH2DX_H +#define _INTERPFROMMESHTOMESH2DX_H + +#include "../../objects/objects.h" + +/* 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_values,int num_default_values,Contour** contours, int numcontours ); + +#endif