Changeset 5032
- Timestamp:
- 08/06/10 14:48:41 (15 years ago)
- Location:
- issm/trunk/src
- Files:
-
- 7 added
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk/src/c/Makefile.am ¶
r4983 r5032 479 479 ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\ 480 480 ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\ 481 ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\ 482 ./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\ 481 483 ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\ 482 484 ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h\ … … 1015 1017 ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\ 1016 1018 ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\ 1019 ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\ 1020 ./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\ 1017 1021 ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\ 1018 1022 ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\ -
TabularUnified issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp ¶
r5016 r5032 8 8 #include "../../toolkits/toolkits.h" 9 9 #include "../../objects/objects.h" 10 #include "../modules.h" 10 11 11 12 using namespace bamg; … … 13 14 14 15 int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data, 15 double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value){16 16 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){ 17 17 18 /*Output*/ 18 19 double* data_interp=NULL; … … 30 31 int verbose=0; 31 32 33 /*default values: */ 34 Vec vec_incontour=NULL; 35 double* incontour=NULL; 36 bool skip_bamg=false; 37 38 32 39 /*Checks*/ 33 40 if (data_cols<=0){ … … 36 43 if (data_rows!=nods_data && data_rows!=nels_data){ 37 44 ISSMERROR("data provided should have either %i or %i lines (not %i)",nods_data,nels_data,data_rows); 45 } 46 if((num_default_values) && (data_cols>1)){ 47 ISSMERROR("data provided can only have 1 column if a default value is provided"); 48 } 49 50 /*If default values supplied, figure out which nodes are inside the contour, including the border of the contour: */ 51 if(num_default_values){ 52 ContourToNodesx( &vec_incontour,x_interp,y_interp,nods_interp,contours,numcontours,1); 53 VecToMPISerial(&incontour,vec_incontour); 38 54 } 39 55 … … 50 66 if (verbose) printf("Loop over the nodes\n"); 51 67 for(i=0;i<nods_interp;i++){ 68 69 /*reset skip_bamg: */ 70 skip_bamg=false; 52 71 53 //Get current point coordinates 54 r.x=x_interp[i]; r.y=y_interp[i]; 55 I2 I=Th.toI2(r); 56 57 //Find triangle holding r/I 58 Triangle &tb=*Th.FindTriangleContaining(I,dete); 59 60 // internal point 61 if (tb.det>0){ 62 //Area coordinate 63 areacoord[0]= (double) dete[0]/ tb.det; 64 areacoord[1]= (double) dete[1] / tb.det; 65 areacoord[2]= (double) dete[2] / tb.det; 66 //3 vertices of the triangle 67 i0=Th.Number(tb[0]); 68 i1=Th.Number(tb[1]); 69 i2=Th.Number(tb[2]); 70 //triangle number 71 it=Th.Number(tb); 72 /*figure out if we should skip bamg logic: */ 73 if(num_default_values){ 74 if(!incontour[i]){ 75 /*This node is not inside the contour. Skip Bamg logic and apply default value.: */ 76 skip_bamg=true; 77 } 72 78 } 73 79 74 //external point 75 else { 76 //Get closest adjacent triangle (inside the mesh) 77 TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj(); 78 int k=ta; 79 Triangle &tc=*(Triangle*)ta; 80 //Area coordinate 81 areacoord[VerticesOfTriangularEdge[k][1]] = aa; 82 areacoord[VerticesOfTriangularEdge[k][0]] = bb; 83 areacoord[OppositeVertex[k]] = 1 - aa -bb; 84 //3 vertices of the triangle 85 i0=Th.Number(tc[0]); 86 i1=Th.Number(tc[1]); 87 i2=Th.Number(tc[2]); 88 //triangle number 89 it=Th.Number(tc); 90 } 80 if(skip_bamg==false){ 91 81 92 /*Last step, P1 interpolation*/ 93 if (data_rows==nods_data){ 94 for (j=0;j<data_cols;j++){ 95 data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j]; 82 //Get current point coordinates 83 r.x=x_interp[i]; r.y=y_interp[i]; 84 I2 I=Th.toI2(r); 85 86 //Find triangle holding r/I 87 Triangle &tb=*Th.FindTriangleContaining(I,dete); 88 89 // internal point 90 if (tb.det>0){ 91 //Area coordinate 92 areacoord[0]= (double) dete[0]/ tb.det; 93 areacoord[1]= (double) dete[1] / tb.det; 94 areacoord[2]= (double) dete[2] / tb.det; 95 //3 vertices of the triangle 96 i0=Th.Number(tb[0]); 97 i1=Th.Number(tb[1]); 98 i2=Th.Number(tb[2]); 99 //triangle number 100 it=Th.Number(tb); 101 } 102 //external point 103 else { 104 //Get closest adjacent triangle (inside the mesh) 105 TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj(); 106 int k=ta; 107 Triangle &tc=*(Triangle*)ta; 108 //Area coordinate 109 areacoord[VerticesOfTriangularEdge[k][1]] = aa; 110 areacoord[VerticesOfTriangularEdge[k][0]] = bb; 111 areacoord[OppositeVertex[k]] = 1 - aa -bb; 112 //3 vertices of the triangle 113 i0=Th.Number(tc[0]); 114 i1=Th.Number(tc[1]); 115 i2=Th.Number(tc[2]); 116 //triangle number 117 it=Th.Number(tc); 118 } 119 120 if (data_rows==nods_data){ 121 for (j=0;j<data_cols;j++){ 122 data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j]; 123 } 124 } 125 else{ 126 for (j=0;j<data_cols;j++){ 127 if (it<0 || it>=nels_data){ 128 ISSMERROR("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data); 129 } 130 data_interp[i*data_cols+j]=data[data_cols*it+j]; 131 } 96 132 } 97 133 } 98 134 else{ 99 for (j=0;j<data_cols;j++){ 100 if (it<0 || it>=nels_data){ 101 ISSMERROR("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data); 102 } 103 data_interp[i*data_cols+j]=data[data_cols*it+j]; 104 } 135 if(num_default_values==1) data_interp[i]=default_values[0]; 136 else data_interp[i]=default_values[i]; 105 137 } 106 107 138 } 108 139 -
TabularUnified issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h ¶
r3913 r5032 10 10 /* local prototypes: */ 11 11 int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data, 12 double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value);12 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 ); 13 13 14 14 #endif -
TabularUnified issm/trunk/src/c/modules/modules.h ¶
r4983 r5032 30 30 #include "./InputDuplicatex/InputDuplicatex.h" 31 31 #include "./InputScalex/InputScalex.h" 32 #include "./InterpFromMesh2dx/InterpFromMesh2dx.h" 32 33 #include "./InterpFromGridToMeshx/InterpFromGridToMeshx.h" 33 34 #include "./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h" -
TabularUnified issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.cpp ¶
r4981 r5032 5 5 6 6 void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ 7 8 int i; 7 9 8 10 /*input: */ … … 26 28 int y_interp_rows; 27 29 28 double default_value; 30 double* default_values=NULL; 31 int num_default_values=0; 32 33 //contours 34 mxArray* matlabstructure=NULL; 35 Contour** contours=NULL; 36 int numcontours; 37 Contour* contouri=NULL; 29 38 30 39 /*Intermediary*/ … … 41 50 42 51 /*checks on arguments on the matlab side: */ 43 CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InterpFromMeshToMesh2dUsage); 52 if(nlhs!=NLHS){ 53 InterpFromMeshToMesh2dUsage(); 54 ISSMERROR("InterpFromMeshToMesh2dUsage usage error"); 55 } 56 if((nrhs!=6) & (nrhs!=8)){ 57 InterpFromMeshToMesh2dUsage(); 58 ISSMERROR("InterpFromMeshToMesh2dUsage usage error"); 59 } 44 60 45 61 /*Input datasets: */ … … 51 67 FetchData(&x_interp,&x_interp_rows,NULL,XINTERPHANDLE); 52 68 FetchData(&y_interp,&y_interp_rows,NULL,YINTERPHANDLE); 53 FetchData(&default_value,DEFAULTHANDLE); 69 70 if(nrhs==8){ 71 72 /*Call expread on filename to build a contour array in the matlab workspace: */ 73 mexCallMATLAB( 1, &matlabstructure, 1, (mxArray**)&FILENAME, "expread"); 74 75 /*default values: */ 76 FetchData(&default_values,&num_default_values,DEFAULTHANDLE); 77 78 /*contours: */ 79 numcontours=mxGetNumberOfElements(matlabstructure); 80 contours=(Contour**)xmalloc(numcontours*sizeof(Contour*)); 81 for(i=0;i<numcontours;i++){ 82 //allocate 83 contouri=(Contour*)xmalloc(sizeof(Contour)); 84 //retrieve dimension of this contour. 85 contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods")); 86 //set pointers. 87 contouri->x=mxGetPr(mxGetField(matlabstructure,i,"x")); 88 contouri->y=mxGetPr(mxGetField(matlabstructure,i,"y")); 89 *(contours+i)=contouri; 90 } 91 92 /* Debugging of contours :{{{1*/ 93 /*for(i=0;i<numcontours;i++){ 94 printf("\nContour echo: contour number %i / %i\n",i+1,numcontours); 95 contouri=*(contours+i); 96 printf(" Number of grids %i\n",contouri->nods); 97 for (j=0;j<contouri->nods;j++){ 98 printf(" %lf %lf\n",*(contouri->x+j),*(contouri->y+j)); 99 } 100 }*/ 101 /*}}}*/ 102 } 103 else{ 104 num_default_values=0; 105 default_values=NULL; 106 numcontours=0; 107 contours=NULL; 108 } 54 109 55 110 /*some checks*/ … … 71 126 /* Run core computations: */ 72 127 if (verbose) printf("Call core\n"); 73 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 );128 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); 74 129 75 130 /*Write data: */ … … 88 143 _printf_("\n"); 89 144 _printf_(" Usage:\n"); 90 _printf_(" data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp,default_value);\n"); 145 _printf_(" data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp);\n"); 146 _printf_(" or data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp,default_value,contourname);\n"); 91 147 _printf_("\n"); 92 148 _printf_(" index: index of the mesh where data is defined\n"); … … 94 150 _printf_(" data: matrix holding the data to be interpolated onto the mesh. (one column per field)\n"); 95 151 _printf_(" x_interp,y_interp: coordinates of the points onto which we interpolate.\n"); 96 _printf_(" default_value: default value if no data is found (holes).\n"); 152 _printf_(" if default_value and contourname not specified: linear interpolation will happen on all x_interp,y_interp.\n"); 153 _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"); 154 _printf_(" note that default_value is either a scalar, or a vector of size length(x_interp)\n"); 97 155 _printf_(" data_interp: vector of mesh interpolated data.\n"); 98 156 _printf_("\n"); -
TabularUnified issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.h ¶
r4236 r5032 25 25 #define YINTERPHANDLE prhs[5] 26 26 #define DEFAULTHANDLE prhs[6] 27 #define FILENAME prhs[7] 27 28 28 29 /* serial output macros: */ … … 32 33 #undef NLHS 33 34 #define NLHS 1 34 #undef NRHS35 #define NRHS 736 35 37 36 #endif -
TabularUnified issm/trunk/src/mex/Makefile.am ¶
r5004 r5032 35 35 InterpFromMeshToMesh3d \ 36 36 InterpFromMeshToGrid \ 37 InterpFromMesh2d \ 37 38 MassFlux\ 38 39 Mergesolutionfromftog\ … … 191 192 InterpFromMeshToGrid/InterpFromMeshToGrid.h 192 193 194 InterpFromMesh2d_SOURCES = InterpFromMesh2d/InterpFromMesh2d.cpp\ 195 InterpFromMesh2d/InterpFromMesh2d.h 196 193 197 AverageFilter_SOURCES = AverageFilter/AverageFilter.cpp\ 194 198 AverageFilter/AverageFilter.h
Note:
See TracChangeset
for help on using the changeset viewer.