Changeset 2362

Timestamp:

10/02/09 10:47:13 (15 years ago)

Author:

Mathieu Morlighem

Message:

Reduced interpolation time for 3d meshes

File:

• issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp (modified) (4 diffs)

issm/trunk/src/c/InterpFromMeshToMesh3dx/InterpFromMeshToMesh3dx.cpp

@@ -22 +22 @@
         double zeta,bed,surface;
         double data_value;
-        double xmin,xmax;
-        double ymin,ymax;
+        double x_prime_min,x_prime_max;
+        double y_prime_min,y_prime_max;
+        double x_tria_min,y_tria_min;
+        double x_tria_max,y_tria_max;
 
         /*some checks*/
@@ -45 +47 @@
 
         /*Get prime mesh extrema coordinates*/
-        xmin=x_prime[0]; xmax=x_prime[0];ymin=y_prime[0]; ymax=y_prime[0];
+        x_prime_min=x_prime[0]; x_prime_max=x_prime[0];y_prime_min=y_prime[0]; y_prime_max=y_prime[0];
         for (i=1;i<nods_prime;i++){
-                if (x_prime[i]<xmin) xmin=x_prime[i];
-                if (x_prime[i]>xmax) xmax=x_prime[i];
-                if (y_prime[i]<ymin) ymin=y_prime[i];
-                if (y_prime[i]>ymax) ymax=y_prime[i];
+                if (x_prime[i]<x_prime_min) x_prime_min=x_prime[i];
+                if (x_prime[i]>x_prime_max) x_prime_max=x_prime[i];
+                if (y_prime[i]<y_prime_min) y_prime_min=y_prime[i];
+                if (y_prime[i]>y_prime_max) y_prime_max=y_prime[i];
         }
 
@@ -64 +66 @@
                 if (debug && fmod((double)i,(double)100)==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)i/nels_data*100);
 
+                /*Get extrema coordinates of current elements*/
+                x_tria_min=x_data[(int)index_data[6*i+0]-1]; x_tria_max=x_tria_min;
+                y_tria_min=y_data[(int)index_data[6*i+0]-1]; y_tria_max=y_tria_min;
+                for (j=1;j<3;j++){
+                        if(x_data[(int)index_data[6*i+j]-1]<x_tria_min) x_tria_min=x_data[(int)index_data[6*i+j]-1];
+                        if(x_data[(int)index_data[6*i+j]-1]>x_tria_max) x_tria_max=x_data[(int)index_data[6*i+j]-1];
+                        if(y_data[(int)index_data[6*i+j]-1]<y_tria_min) y_tria_min=y_data[(int)index_data[6*i+j]-1];
+                        if(y_data[(int)index_data[6*i+j]-1]>y_tria_max) y_tria_max=y_data[(int)index_data[6*i+j]-1];
+                }
+
                 /*if there is no point inside the domain, go to next iteration*/
-                if ( (x_data[(int)index_data[6*i+0]-1]<xmin) && (x_data[(int)index_data[6*i+1]-1]<xmin) && (x_data[(int)index_data[6*i+2]-1]<xmin)) continue;
-                if ( (x_data[(int)index_data[6*i+0]-1]>xmax) && (x_data[(int)index_data[6*i+1]-1]>xmax) && (x_data[(int)index_data[6*i+2]-1]>xmax)) continue;
-                if ( (y_data[(int)index_data[6*i+0]-1]<ymin) && (y_data[(int)index_data[6*i+1]-1]<ymin) && (y_data[(int)index_data[6*i+2]-1]<ymin)) continue;
-                if ( (y_data[(int)index_data[6*i+0]-1]>ymax) && (y_data[(int)index_data[6*i+1]-1]>ymax) && (y_data[(int)index_data[6*i+2]-1]>ymax)) continue;
+                if ( x_prime_max < x_tria_min ) continue; 
+                if ( x_prime_min > x_tria_max ) continue; 
+                if ( y_prime_max < y_tria_min ) continue; 
+                if ( y_prime_min > y_tria_max ) continue; 
 
                 /*get area of the current element (Jacobian = 2 * area)*/
@@ -78 +90 @@
                 /*loop over the prime nodes*/
                 for (j=0;j<nods_prime;j++){
+
+                        /*if the current point is not in the triangle, continue*/
+                        if ( x_prime[j] < x_tria_min ) continue; 
+                        if ( x_prime[j] > x_tria_max ) continue; 
+                        if ( y_prime[j] < y_tria_min ) continue; 
+                        if ( y_prime[j] > y_tria_max ) continue; 
 
                         /*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/