Changeset 12604

Timestamp:

07/03/12 09:17:57 (13 years ago)

Author:

Mathieu Morlighem

Message:

only interpolate for points that are inside the model domain

File:

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

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

@@ -51 +51 @@
         /*read background mesh*/
         Mesh Th(index_data,x_data,y_data,nods_data,nels_data); 
+
+        /*Get reference number (for subdomains)*/
+        long* reft = xNew<long>(Th.nbt);
+        Th.TriangleReferenceList(reft);
         Th.CreateSingleVertexToTriangleConnectivity();
 
@@ -56 +60 @@
         for(i=0;i<N_interp;i++){
 
-                //Get current point coordinates
+                /*Get current point coordinates*/
                 r.x=x_interp[i]; r.y=y_interp[i];
                 I2 I=Th.R2ToI2(r);
 
-                //Find triangle holding r/I
+                /*Find triangle holding r/I*/
                 Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
 
-                // internal point 
+                /*point inside convex*/
                 if (tb.det>0){ 
-                        //Area coordinate
+
+                        /*Area coordinates*/
                         areacoord[0]= (double) dete[0]/tb.det;
                         areacoord[1]= (double) dete[1]/tb.det;
                         areacoord[2]= (double) dete[2]/tb.det;
-                        //3 vertices of the triangle
+                        /*3 vertices of the triangle*/
                         i0=Th.GetId(tb[0]);
                         i1=Th.GetId(tb[1]);
                         i2=Th.GetId(tb[2]);
-                        //triangle number
+                        /*triangle number*/
                         it=Th.GetId(tb);
+
+                        /*Inside convex but outside mesh*/
+                        if (reft[it]<0 & isdefault){
+                                for(j=0;j<N_data;j++) data_interp[i*N_data+j]=defaultvalue;
+                                continue;
+                        }
                 }
                 //external point
@@ -115 +126 @@
         }
 
-        /*Assign output pointers:*/
+        /*clean-up and return*/
+        xDelete<long>(reft);
         *pdata_interp=data_interp;
-
-        /*No error return*/
         return 1;
 }