Changeset 23085

Timestamp:

08/10/18 16:30:54 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: fixing an indexing problem and sped up distance test

File:

• issm/trunk-jpl/src/c/classes/FemModel.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1050 +1050 @@
         /*3: Add distance input to all elements*/
         IssmDouble* distances = xNew<IssmDouble>(vertices->Size());
-        IssmDouble  d,xn,yn;
-        for(int i=0;i<vertices->Size();i++){
-                Vertex* vertex=dynamic_cast<Vertex*>(this->vertices->GetObjectByOffset(i));
-                IssmDouble dmin = 1.e+50;
+        IssmDouble  d,xn,yn,dmin;
+        int         last = -1;
+        for(int v=0;v<vertices->Size();v++){
+
+                Vertex* vertex=dynamic_cast<Vertex*>(this->vertices->GetObjectByOffset(v));
+                IssmDouble x = vertex->x;
+                IssmDouble y = vertex->y;
+
+                /*Most of the time the last checked segment is the closest so start with that one*/
+                if(last>0){
+                        dmin = sqrt(pow(allsegmentlist[4*last+0] - x,2) + pow(y-allsegmentlist[4*last+1],2));
+                }
+                else{
+                        dmin = 1.e+50;
+                }
 
                 for(int i=0;i<numseg;i++){
-                        IssmDouble x = vertex->x;
-                        IssmDouble y = vertex->y;
 
                         /*Skip if tip is more than 10xdmin away*/
-                        if( pow(allsegmentlist[4*i+0] - x,2) + pow(y-allsegmentlist[4*i+1],2) > pow(10*dmin,2)) continue;
+                        if(fabs(allsegmentlist[4*i+0] - x)>10*dmin) continue;
+                        if(fabs(allsegmentlist[4*i+0] - y)>10*dmin) continue;
 
                         IssmDouble l2 = (allsegmentlist[4*i+2]-allsegmentlist[4*i+0])*(allsegmentlist[4*i+2]-allsegmentlist[4*i+0]) + (allsegmentlist[4*i+3]-allsegmentlist[4*i+1])*(allsegmentlist[4*i+3]-allsegmentlist[4*i+1]);
@@ -1067 +1077 @@
                         if(l2==0.){
                                 d = (x-allsegmentlist[4*i+0])*(x-allsegmentlist[4*i+0])+(y-allsegmentlist[4*i+1])*(y-allsegmentlist[4*i+1]);
-                                if(d<dmin) dmin = d;
+                                if(d<dmin){
+                                        dmin = d;
+                                        last = i;
+                                }
                                 continue;
                         }
@@ -1090 +1103 @@
                         }
 
-                        if(d<dmin) dmin = d;
+                        if(d<dmin){
+                                dmin = d;
+                                last = i;
+                        }
                 }