Changeset 18788

Timestamp:

11/14/14 16:42:44 (10 years ago)

Author:

Eric.Larour

Message:

CHG: improvements to this module, this could seriously impact some of the nightly runs dealing with UQ, so to be checked.

Location:

issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx

Files:

• MeshProfileIntersectionx.cpp (modified) (5 diffs)
• MeshProfileIntersectionx.h (modified) (1 diff)

issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp ¶

@@ -129 +129 @@
                 ysegment[0]=yc[i];
                 ysegment[1]=yc[i+1];
-                ElementSegment(segments_dataset,el, xnodes,ynodes,xsegment,ysegment);
-        }
-}/*}}}*/
-void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment){/*{{{*/
+                /*if (el==318 && i==9){
+                        _printf_("contour: " << i << " " << xsegment[0] << " " << ysegment[0] << " " << xsegment[1] << " " << ysegment[1] 
+                                << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << 
+                                ynodes[2] << "\n");
+                }*/
+                ElementSegment(segments_dataset,el, i, xnodes,ynodes,xsegment,ysegment);
+        }
+}/*}}}*/
+void ElementSegment(DataSet* segments_dataset,int el, int contouri, double* xnodes,double* ynodes,double* xsegment,double* ysegment){/*{{{*/
 
         /*We have a tria element (xnodes,ynodes) and a segment (xsegment,ysegment). Find whether they intersect. 
@@ -159 +164 @@
         edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
 
-        /*edge can be either IntersectEnum (one iand only one intersection between the edge and the segment), ColinearEnum (edge and segment are collinear) and SeparateEnum (no intersection): */
+        /*edge can be either IntersectEnum (one and only one intersection between the edge and the segment), ColinearEnum (edge and segment are collinear) and SeparateEnum (no intersection): */
+                
+        /*if (el==318 && contouri==9){
+                _printf_(edge1 << " " << edge2 << " " << edge3 << " "  << alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " "  << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
+        
+        _printf_("Bool" << (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum));
+        }*/
 
         if(    (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum)   ){
-                /*This case is impossible: */
-                _error_("error: a line cannot go through 3 different vertices!");
+
+                /*This can only be the case if the segment intersected through one vertex, meaning a pair from alpha1, beta1 or gamma1  is 0:*/
+                if (alpha1!=0 && alpha1!=1){
+                        /*The vertex opposite edge 1 was intersected:*/
+                        xfinal[0]=xsegment[0]+alpha1*(xsegment[1]-xsegment[0]);
+                        yfinal[0]=ysegment[0]+alpha1*(ysegment[1]-ysegment[0]);
+                        xfinal[1]=xnodes[2];
+                        yfinal[1]=ynodes[2];
+                }
+                else if (beta1!=0 && beta1!=1){
+                        /*The vertex opposite edge 2 was intersected:*/
+                        xfinal[0]=xsegment[0]+beta1*(xsegment[1]-xsegment[0]);
+                        yfinal[0]=ysegment[0]+beta1*(ysegment[1]-ysegment[0]);
+                        xfinal[1]=xnodes[0];
+                        yfinal[1]=ynodes[0];
+                }
+                else if (gamma1!=0 && gamma1!=1){
+                        /*The vertex opposite edge 3 was intersected:*/
+                        xfinal[0]=xsegment[0]+gamma1*(xsegment[1]-xsegment[0]);
+                        yfinal[0]=ysegment[0]+gamma1*(ysegment[1]-ysegment[0]);
+                        xfinal[1]=xnodes[1];
+                        yfinal[1]=ynodes[1];
+                }
+                segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
+
+                
+                /*This case is impossible: not quite! */
+                //_printf_(alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " "  << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
+                /* _error_("error: a line cannot go through 3 different vertices!");*/
         }
         else if(    ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum))   ){
@@ -187 +225 @@
         }
         else if(  (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)   ){
+        
+                /*if (el==318 && contouri==9){
+                        _printf_("hello" <<  " NodeInElement 0 " << (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) <<  " NodeInElement 1 " << (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])));
+                }*/
 
                 /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element, 
@@ -196 +238 @@
                         if(edge3==IntersectEnum){coord2=gamma1;}
                 }
-                else{
+                else if (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])){
                         if(edge1==IntersectEnum){coord1=alpha1;}
                         if(edge2==IntersectEnum){coord1=beta1;}
                         if(edge3==IntersectEnum){coord1=gamma1;}
                         coord2=1.0;
+                }
+                else{
+                        double tolerance=1e-10;
+                        /*Ok, we have an issue here. Probably one of the segments' end is on a vertex, within a certain tolerance!*/
+                        if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[0],ysegment[0],tolerance) ||
+                                IsIdenticalNode(xnodes[1],ynodes[1],xsegment[0],ysegment[0],tolerance) ||
+                                IsIdenticalNode(xnodes[2],ynodes[2],xsegment[0],ysegment[0],tolerance)){
+                                
+                                /*ok, segments[0] is common to one of our vertices: */
+                                //if (el==318 && contouri==9){ _printf_("ok1" << "\n"); }
+                                coord1=0;
+                                if(edge1==IntersectEnum){coord2=alpha1;}
+                                if(edge2==IntersectEnum){coord2=beta1;}
+                                if(edge3==IntersectEnum){coord2=gamma1;}
+                        }
+                        else if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[1],ysegment[1],tolerance) ||
+                                     IsIdenticalNode(xnodes[1],ynodes[1],xsegment[1],ysegment[1],tolerance) ||
+                                     IsIdenticalNode(xnodes[2],ynodes[2],xsegment[1],ysegment[1],tolerance)){
+
+                                /*ok, segments[1] is common to one of our vertices: */
+                                //if (el==318 && contouri==9){ _printf_("ok2" << "\n"); }
+                                if(edge1==IntersectEnum){coord1=alpha1;}
+                                if(edge2==IntersectEnum){coord1=beta1;}
+                                if(edge3==IntersectEnum){coord1=gamma1;}
+                                coord2=1.0;
+                        }
                 }
 
@@ -344 +412 @@
         return IntersectEnum;
 }/*}}}*/
+bool IsIdenticalNode(double x1, double y1, double x2, double y2, double tolerance){ /*{{{*/
+
+        if (sqrt(pow(x1-x2,2.0) + pow(y1-y2,2))<tolerance)return true;
+        else return false;
+
+}/*}}}*/

issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h ¶

@@ -13 +13 @@
 void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes);
 void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes);
-void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment);
+void ElementSegment(DataSet* segments_dataset,int el,int contouri, double* xnodes,double* ynodes,double* xsegment,double* ysegment);
 int  SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2);
 bool NodeInElement(double* xnodes, double* ynodes, double x, double y);
+bool IsIdenticalNode(double x1, double y1, double x2, double y2, double tolerance);
 
 #endif /* _MESHPROFILEINTERSECTIONX_H */