Changeset 22729

Timestamp:

05/01/18 02:28:49 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: convert index to int

Location:

issm/trunk-jpl/src

Files:

• c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp (modified) (7 diffs)
• c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.h (modified) (1 diff)
• wrappers/InterpFromMeshToGrid/InterpFromMeshToGrid.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp ¶

@@ -6 +6 @@
 #include "../../shared/shared.h"
 
-void InterpFromMeshToGridx(double** pgriddata,double* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh,int data_length,double* x_grid,double* y_grid,int nlines,int ncols,double default_value){
+void InterpFromMeshToGridx(double** pgriddata,int* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh,int data_length,double* x_grid,double* y_grid,int nlines,int ncols,double default_value){
 
         /*Intermediary*/
-        int    i,j;
         int    i1,i2,j1,j2;
         int    interpolation_type;
@@ -28 +27 @@
         else _error_("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
 
-        /*First, allocate pointers: */
-        double* griddata=xNewZeroInit<double>(nlines*ncols);
+        /*First, prepare output*/
+        double* griddata=xNew<double>(nlines*ncols);
+        for(int i=0;i<nlines;i++) for(int j=0;j<ncols; j++) griddata[i*ncols+j]=default_value;
 
-        /*Set debug to 1 if there are lots of elements*/
+        /*Set debug to "true" if there are lots of elements*/
         bool debug=(bool)((double)ncols*nlines*nels >= 5.e10);
-
-        /*Initialize coordintes and griddata*/
-        for(i=0;i<nlines;i++) for(j=0;j<ncols; j++) griddata[i*ncols+j]=default_value;
 
         /*figure out if x or y are flipped*/
         bool xflip = false;
         bool yflip = false;
-        if (x_grid[1]-x_grid[0]<0) xflip=true;
-        if (y_grid[1]-y_grid[0]<0) yflip=true;
+        if(x_grid[1]-x_grid[0]<0) xflip=true;
+        if(y_grid[1]-y_grid[0]<0) yflip=true;
 
         /*Get min/max coordinates of the grid*/
@@ -67 +64 @@
 
                 /*display current iteration*/
-                if(debug && fmod((double)n,(double)100)==0)
+                if(debug && fmod(double(n),1000)==0)
                  _printf_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(n)/double(nels)*100<<"%   ");
 
                 /*Get extrema coordinates of current elements*/
-                x_tria_min=x_mesh[(int)index_mesh[3*n+0]-1]; x_tria_max=x_tria_min;
-                y_tria_min=y_mesh[(int)index_mesh[3*n+0]-1]; y_tria_max=y_tria_min;
-                for(i=1;i<3;i++){
-                        if(x_mesh[(int)index_mesh[3*n+i]-1]<x_tria_min) x_tria_min=x_mesh[(int)index_mesh[3*n+i]-1];
-                        if(x_mesh[(int)index_mesh[3*n+i]-1]>x_tria_max) x_tria_max=x_mesh[(int)index_mesh[3*n+i]-1];
-                        if(y_mesh[(int)index_mesh[3*n+i]-1]<y_tria_min) y_tria_min=y_mesh[(int)index_mesh[3*n+i]-1];
-                        if(y_mesh[(int)index_mesh[3*n+i]-1]>y_tria_max) y_tria_max=y_mesh[(int)index_mesh[3*n+i]-1];
+                x_tria_min=x_mesh[index_mesh[3*n+0]-1]; x_tria_max=x_tria_min;
+                y_tria_min=y_mesh[index_mesh[3*n+0]-1]; y_tria_max=y_tria_min;
+                for(int i=1;i<3;i++){
+                        if(x_mesh[index_mesh[3*n+i]-1]<x_tria_min) x_tria_min=x_mesh[index_mesh[3*n+i]-1];
+                        if(x_mesh[index_mesh[3*n+i]-1]>x_tria_max) x_tria_max=x_mesh[index_mesh[3*n+i]-1];
+                        if(y_mesh[index_mesh[3*n+i]-1]<y_tria_min) y_tria_min=y_mesh[index_mesh[3*n+i]-1];
+                        if(y_mesh[index_mesh[3*n+i]-1]>y_tria_max) y_tria_max=y_mesh[index_mesh[3*n+i]-1];
                 }
 
@@ -83 +80 @@
                 if( (x_tria_min>x_grid_max) || (x_tria_max<x_grid_min) || (y_tria_min>y_grid_max) || (y_tria_max<y_grid_min) ) continue;
 
-                /*Get indices i and j that form a square around the currant triangle*/
+                /*Get indices i and j that form a square around the current triangle*/
                 if(yflip){
                         i1=max(0,       (int)floor((y_tria_max-y_grid_max)/yposting)-1);
@@ -103 +100 @@
                 /*get area of the current element (Jacobian = 2 * area)*/
                 //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
-                area=x_mesh[(int)index_mesh[3*n+1]-1]*y_mesh[(int)index_mesh[3*n+2]-1]-y_mesh[(int)index_mesh[3*n+1]-1]*x_mesh[(int)index_mesh[3*n+2]-1]
-                  +  x_mesh[(int)index_mesh[3*n+0]-1]*y_mesh[(int)index_mesh[3*n+1]-1]-y_mesh[(int)index_mesh[3*n+0]-1]*x_mesh[(int)index_mesh[3*n+1]-1]
-                  +  x_mesh[(int)index_mesh[3*n+2]-1]*y_mesh[(int)index_mesh[3*n+0]-1]-y_mesh[(int)index_mesh[3*n+2]-1]*x_mesh[(int)index_mesh[3*n+0]-1];
+                area=x_mesh[index_mesh[3*n+1]-1]*y_mesh[index_mesh[3*n+2]-1]-y_mesh[index_mesh[3*n+1]-1]*x_mesh[index_mesh[3*n+2]-1]
+                  +  x_mesh[index_mesh[3*n+0]-1]*y_mesh[index_mesh[3*n+1]-1]-y_mesh[index_mesh[3*n+0]-1]*x_mesh[index_mesh[3*n+1]-1]
+                  +  x_mesh[index_mesh[3*n+2]-1]*y_mesh[index_mesh[3*n+0]-1]-y_mesh[index_mesh[3*n+2]-1]*x_mesh[index_mesh[3*n+0]-1];
 
                 /*Go through x_grid and y_grid and interpolate if necessary*/
-                for(i=i1;i<=i2;i++){
+                for(int i=i1;i<=i2;i++){
 
                         //exit if y not between y_tria_min and y_tria_max
                         if((y_grid[i]>y_tria_max) || (y_grid[i]<y_tria_min)) continue;
 
-                        for(j=j1;j<=j2; j++){
+                        for(int j=j1;j<=j2; j++){
 
                                 //exit if x not between x_tria_min and x_tria_max
@@ -119 +116 @@
 
                                 /*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
-                                area_1=((x_grid[j]-x_mesh[(int)index_mesh[3*n+2]-1])*(y_mesh[(int)index_mesh[3*n+1]-1]-y_mesh[(int)index_mesh[3*n+2]-1]) 
-                                                        -  (y_grid[i]-y_mesh[(int)index_mesh[3*n+2]-1])*(x_mesh[(int)index_mesh[3*n+1]-1]-x_mesh[(int)index_mesh[3*n+2]-1]))/area;
+                                area_1=((x_grid[j]-x_mesh[index_mesh[3*n+2]-1])*(y_mesh[index_mesh[3*n+1]-1]-y_mesh[index_mesh[3*n+2]-1]) 
+                                                        -  (y_grid[i]-y_mesh[index_mesh[3*n+2]-1])*(x_mesh[index_mesh[3*n+1]-1]-x_mesh[index_mesh[3*n+2]-1]))/area;
                                 /*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
-                                area_2=((x_mesh[(int)index_mesh[3*n+0]-1]-x_mesh[(int)index_mesh[3*n+2]-1])*(y_grid[i]-y_mesh[(int)index_mesh[3*n+2]-1]) 
-                                                        - (y_mesh[(int)index_mesh[3*n+0]-1]-y_mesh[(int)index_mesh[3*n+2]-1])*(x_grid[j]-x_mesh[(int)index_mesh[3*n+2]-1]))/area;
+                                area_2=((x_mesh[index_mesh[3*n+0]-1]-x_mesh[index_mesh[3*n+2]-1])*(y_grid[i]-y_mesh[index_mesh[3*n+2]-1]) 
+                                                        - (y_mesh[index_mesh[3*n+0]-1]-y_mesh[index_mesh[3*n+2]-1])*(x_grid[j]-x_mesh[index_mesh[3*n+2]-1]))/area;
                                 /*Get third area coordinate = 1-area1-area2*/
                                 area_3=1.-area_1-area_2;
 
                                 /*is the current point in the current element?*/
-                                if(area_1>-10e-12 && area_2>-10e-12 && area_3>-10e-12){
+                                if(area_1>-1e-12 && area_2>-1e-12 && area_3>-1e-12){
 
                                         /*Yes ! compute the value on the point*/
                                         if(interpolation_type==1){
                                                 /*nodal interpolation*/
-                                                data_value=area_1*data_mesh[(int)index_mesh[3*n+0]-1]+area_2*data_mesh[(int)index_mesh[3*n+1]-1]+area_3*data_mesh[(int)index_mesh[3*n+2]-1];
+                                                data_value=area_1*data_mesh[index_mesh[3*n+0]-1]+area_2*data_mesh[index_mesh[3*n+1]-1]+area_3*data_mesh[index_mesh[3*n+2]-1];
                                         }
                                         else{
@@ -147 +144 @@
                 }
         }
-        if(debug) _printf_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%  \n");
+        if(debug) _printf_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<fixed<<100.<<"%   \n");
 
         /*Assign output pointers:*/

issm/trunk-jpl/src/c/modules/InterpFromMeshToGridx/InterpFromMeshToGridx.h ¶

@@ -8 +8 @@
 #include "../../toolkits/toolkits.h"
 
-void InterpFromMeshToGridx(double** pgriddata,double* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh,int data_length,double* x_grid,double* y_grid,int nlines,int ncols,double default_value);
+void InterpFromMeshToGridx(double** pgriddata,int* index_mesh, double* x_mesh, double* y_mesh, int nods,int nels, double* data_mesh,int data_length,double* x_grid,double* y_grid,int nlines,int ncols,double default_value);
 
 #endif /* _INTERPFROMMESHTOGRIDX_H */

issm/trunk-jpl/src/wrappers/InterpFromMeshToGrid/InterpFromMeshToGrid.cpp ¶

@@ -23 +23 @@
 
         /*inputs */
-        double* index=NULL;
+        int*    index=NULL;
         double* x=NULL;
         double* y=NULL;
@@ -63 +63 @@
 
         /*Free ressources: */
-        xDelete<double>(index);
+        xDelete<int>(index);
         xDelete<double>(x);
         xDelete<double>(y);