Changeset 21802 for issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp

Timestamp:

07/17/17 21:18:32 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on bamg amr

File:

• issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp

@@ -24 +24 @@
         R2     r;
         I2     I;
-        int    i,j;
+        int    i,j,k;
         int    it;
         int    i0,i1,i2;
@@ -64 +64 @@
                 if(y_data[i]<ymin) ymin=y_data[i];
                 if(y_data[i]>ymax) ymax=y_data[i];
+        }
+
+        /*Create Single vertex to element connectivity*/
+        int* connectivity = xNew<int>(nods_data);
+        for(i=0;i<nels_data;i++){
+                for(j=0;j<3;j++){
+                        k = index_data[i*3+j]-1;
+                        _assert_(k>=0 & k<nods_data);
+                        connectivity[k]=i;
+                }
         }
 
@@ -137 +147 @@
                         /*If we fall outside of the convex or outside of the mesh, return NaN*/
                         if(tb.det<0 || reft[it]<0){
-                                for (j=0;j<N_data;j++){
-                                        data_interp[i*N_data+j]=NAN;
-                                }
+                                _assert_(i0>=0 & i0<nods_data);
+                                it=connectivity[i0]; //or i1 or i2
+                                _assert_(it>=0 && it<nels_data);
+                                for(j=0;j<N_data;j++) data_interp[i*N_data+j]=data[N_data*it+j];
                         }
                         else{
+                                /*Inside the mesh!*/
                                 if(it<0 || it>=nels_data){
                                         _error_("Triangle number " << it << " not in [0 " << nels_data
                                                                 << "], report bug to developers (interpolation point: " <<x_interp[i]<<" "<<y_interp[i]<<")");
                                 }
-                                for (j=0;j<N_data;j++){
-
-                                        data_interp[i*N_data+j]=data[N_data*it+j];
-                                }
+                                for (j=0;j<N_data;j++) data_interp[i*N_data+j]=data[N_data*it+j];
                         }
                 }
@@ -158 +167 @@
         delete Th;
         xDelete<long>(reft);
+        xDelete<int>(connectivity);
         *pdata_interp=data_interp;
         return 1;