Index: /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 18787) +++ /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 18788) @@ -129,8 +129,13 @@ 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,9 +164,42 @@ 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(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,4 +225,8 @@ } 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,9 +238,35 @@ 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,2 +412,8 @@ 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))