Changeset 8157

Timestamp:

05/06/11 08:16:05 (14 years ago)

Author:

Mathieu Morlighem

Message:

Fix problem if matrix is smaller than model mesh

File:

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

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

@@ -27 +27 @@
         double x_grid_max,y_grid_max;
         double data_value;
-        double cornersouth,cornerwest;
+        double cornersouth;
         double* x_grid=NULL;
         double* y_grid=NULL;
@@ -66 +66 @@
         // EAST = ACTUAL WEST !!!!!!!!!!!!!!!
         cornersouth=cornernorth-nlines*yposting;
-        cornerwest =cornereast +ncols *xposting;
         for(i=0;i<nlines;i++) y_grid[i]= cornersouth + yposting/2 + yposting*i;
         for(i=0;i<ncols; i++) x_grid[i]= cornereast  + xposting/2 + xposting*i;
-        for(i=0;i<=nlines;i++)   y_m[i]   = cornersouth + yposting*i;
-        for(i=0;i<=ncols; i++)   x_m[i]   = cornereast  + xposting*i;
+        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*/
@@ -115 +114 @@
                 /*Get indices i and j that form a square around the currant triangle*/
                 if (yflip){
-                        i1=(int)floor((y_tria_max-y_grid_max)/yposting)-1;
-                        i2= (int)ceil((y_tria_min-y_grid_max)/yposting);
+                        i1=max(0,      (int)floor((y_tria_max-y_grid_max)/yposting)-1);
+                        i2=min(nlines, (int)ceil((y_tria_min-y_grid_max)/yposting));
                 }
                 else{
-                        i1=(int)floor((y_tria_min-y_grid_min)/yposting)-1;
-                        i2= (int)ceil((y_tria_max-y_grid_min)/yposting);
+                        i1=max(0,      (int)floor((y_tria_min-y_grid_min)/yposting)-1);
+                        i2=min(nlines, (int)ceil((y_tria_max-y_grid_min)/yposting));
                 }
                 if (xflip){
-                        j1=(int)floor((x_tria_max-x_grid_max)/xposting)-1;
-                        j2= (int)ceil((x_tria_min-x_grid_max)/xposting);
+                        j1=max(0,     (int)floor((x_tria_max-x_grid_max)/xposting)-1);
+                        j2=min(ncols, (int)ceil((x_tria_min-x_grid_max)/xposting));
                 }
                 else{
-                        j1=(int)floor((x_tria_min-x_grid_min)/xposting)-1;
-                        j2= (int)ceil((x_tria_max-x_grid_min)/xposting);
+                        j1=max(0,     (int)floor((x_tria_min-x_grid_min)/xposting)-1);
+                        j2=min(ncols, (int)ceil((x_tria_max-x_grid_min)/xposting));
                 }
-                /*Check indices*/
-                if(i1<0)      _error_("first index is 0!");
-                if(j1<0)      _error_("first index is 0!");
-                if(i2>=nlines)_error_("index exceeds matrix dimension! (%i >= nlines  %i)",i2,nlines);
-                if(j2>=ncols) _error_("index exceeds matrix dimension! (%i >= ncols   %i)",j2,ncols);
 
                 /*get area of the current element (Jacobian = 2 * area)*/