Changeset 2304

Timestamp:

09/24/09 09:16:11 (15 years ago)

Author:

Mathieu Morlighem

Message:

accelerated interpolation process

File:

• issm/trunk/src/c/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp (modified) (6 diffs)

issm/trunk/src/c/InterpFromMeshToGridx/InterpFromMeshToGridx.cpp

@@ -18 +18 @@
         /*Intermediary*/
         int    i,j,n;
+        int    i1,i2,j1,j2;
         int    interpolation_type;
         bool   debug;
+        int    xflip,yflip;
         double area;
         double area_1,area_2,area_3;
@@ -57 +59 @@
         debug=(bool)((double)ncols*nlines*nels >= pow((double)10,(double)10));
 
-        /*Initialize coordintes and griddata*/
+        /*figure out if x or y are flipped*/
+        if (xposting<0) xflip=1;
+        else xflip=0;
+        if (yposting<0) yflip=1;
+        else yflip=0;
+
+        /*Compute coordinates lists*/
         // EAST = ACTUAL WEST !!!!!!!!!!!!!!!
         cornersouth=cornernorth-nlines*yposting;
@@ -65 +73 @@
         for(i=0;i<=nlines;i++)   y_m[i]   = cornersouth + yposting*i;
         for(i=0;i<=ncols; i++)   x_m[i]   = cornereast  + xposting*i;
+
+        /*Initialize coordintes and griddata*/
         for(i=0;i<nlines;i++){
                 for(j=0;j<ncols; j++){
@@ -72 +82 @@
         
         /*Get extreme coordinates of the grid*/
-        if (xposting>0){
-                x_grid_min=cornereast;
-                x_grid_max=cornerwest;
+        if (xflip){
+                x_grid_min=x_grid[ncols-1]; x_grid_max=x_grid[0];
         }
         else{
-                x_grid_min=cornerwest;
-                x_grid_max=cornereast;
+                x_grid_min=x_grid[0]; x_grid_max=x_grid[ncols-1];
         }
-        if (yposting>0){
-                y_grid_min=cornersouth;
-                y_grid_max=cornernorth;
+        if (yflip){
+                y_grid_min=y_grid[nlines-1]; y_grid_max=y_grid[0];
         }
         else{
-                y_grid_min=cornernorth;
-                y_grid_max=cornersouth;
+                y_grid_min=y_grid[0]; y_grid_max=y_grid[nlines-1];
         }
 
@@ -109 +115 @@
                 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*/
+                if (yflip){
+                        i1=floor((y_tria_max-y_grid_max)/yposting)-1;
+                        i2= ceil((y_tria_min-y_grid_max)/yposting);
+                }
+                else{
+                        i1=floor((y_tria_min-y_grid_min)/yposting)-1;
+                        i2= ceil((y_tria_max-y_grid_min)/yposting);
+                }
+                if (xflip){
+                        j1=floor((x_tria_max-x_grid_max)/xposting)-1;
+                        j2= ceil((x_tria_min-x_grid_max)/xposting);
+                }
+                else{
+                        j1=floor((x_tria_min-x_grid_min)/xposting)-1;
+                        j2= ceil((x_tria_max-x_grid_min)/xposting);
+                }
+
                 /*get area of the current element (Jacobian = 2 * area)*/
                 //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
@@ -115 +139 @@
                   +  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];
 
-
                 /*Go through x_grid and y_grid and interpolate if necessary*/
-                for (i=0;i<nlines;i++){
+                for (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=0;j<ncols; j++){
+                        for(j=j1;j<=j2; j++){
 
                                 //exit if x not between x_tria_min and x_tria_max