Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 14217) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 14218) @@ -781,9 +781,4 @@ ./modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp\ ./modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h\ - ./modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp\ - ./modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp\ - ./modules/MeshProfileIntersectionx/ElementSegment.cpp\ - ./modules/MeshProfileIntersectionx/SegmentIntersect.cpp\ - ./modules/MeshProfileIntersectionx/NodeInElement.cpp\ ./modules/ContourToMeshx/ContourToMeshx.cpp\ ./modules/ContourToMeshx/ContourToMeshxt.cpp\ Index: sm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp (revision 14217) +++ (revision ) @@ -1,90 +1,0 @@ -/*! \file ElementSegment.cpp - */ - -#include "./MeshProfileIntersectionx.h" - -void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment){ - - /*We have a tria element (xnodes,ynodes) and a segment (xsegment,ysegment). Find whether they intersect. - * If they do, create a Segment object with the intersection, and add to segments_dataset dataset: */ - - double alpha1,alpha2; - double beta1,beta2; - double gamma1,gamma2; - - int edge1,edge2,edge3; - - double xel[2],yel[2]; - double coord1,coord2; - double xfinal[2],yfinal[2]; - - /*edge 1: */ - xel[0]=xnodes[0]; yel[0]=ynodes[0]; xel[1]=xnodes[1]; yel[1]=ynodes[1]; - edge1=SegmentIntersect(&alpha1,&alpha2, xel,yel,xsegment,ysegment); //alpha1: segment coordinate of intersection. alpha2: same thing for second interesection if it exists (colinear edges) - - /*edge 2: */ - xel[0]=xnodes[1]; yel[0]=ynodes[1]; xel[1]=xnodes[2]; yel[1]=ynodes[2]; - edge2=SegmentIntersect(&beta1,&beta2, xel,yel,xsegment,ysegment); - - /*edge 3: */ - xel[0]=xnodes[2]; yel[0]=ynodes[2]; xel[1]=xnodes[0]; yel[1]=ynodes[0]; - 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): */ - - if( (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum) ){ - /*This case is impossible: */ - _error_("error: a line cannot go through 3 different vertices!"); - } - else if( ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum)) ){ - - /*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */ - if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);} - if((edge2==IntersectEnum) && (edge3==IntersectEnum)) {coord1=min(beta1,gamma1); coord2=max(beta1,gamma1);} - if((edge3==IntersectEnum) && (edge1==IntersectEnum)) {coord1=min(gamma1,alpha1); coord2=max(gamma1,alpha1);} - - /*check this segment did not intersect at a vertex of the tria: */ - if(coord1!=coord2){ - - xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]); - xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]); - yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]); - yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]); - - segments_dataset->AddObject(new Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1])); - } - else{ - /*the segment intersected at the vertex, do not bother with this "0" length segment!:*/ - } - } - else if( (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum) ){ - - /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element, - * this will decide the coordinate: */ - if (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])){ - coord1=0; - if(edge1==IntersectEnum){coord2=alpha1;} - if(edge2==IntersectEnum){coord2=beta1;} - if(edge3==IntersectEnum){coord2=gamma1;} - } - else{ - if(edge1==IntersectEnum){coord1=alpha1;} - if(edge2==IntersectEnum){coord1=beta1;} - if(edge3==IntersectEnum){coord1=gamma1;} - coord2=1.0; - } - - xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]); - xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]); - yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]); - yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]); - - segments_dataset->AddObject(new Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1])); - } - else{ - /*No interesections, but the segment might be entirely inside this triangle!: */ - if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){ - segments_dataset->AddObject(new Segment(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1])); - } - } -} Index: sm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp (revision 14217) +++ (revision ) @@ -1,23 +1,0 @@ -/*! \file ElementSegmentsIntersection.cpp - */ - -#include "./MeshProfileIntersectionx.h" - -void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes){ - - int i; - double xsegment[2]; - double ysegment[2]; - - /*Loop through contour: */ - for(i=0;i** contours,int numcontours){ +void MeshProfileIntersectionx(int** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods, Contour** contours,int numcontours){/*{{{*/ int i,j,k; @@ -15,16 +15,16 @@ /*output: */ - double* segments=NULL; - int numsegs; + int* segments=NULL; + int numsegs; /*intermediary: */ - double** allsegments=NULL; - double* segmentsi=NULL; - int* allnumsegs=NULL; - int numsegsi; - int count; + int** allsegments=NULL; + int* segmentsi=NULL; + int* allnumsegs=NULL; + int numsegsi; + int count; /*Allocate: */ - allsegments=xNew(numcontours); + allsegments=xNew(numcontours); allnumsegs=xNew(numcontours); @@ -51,5 +51,5 @@ /*Out of all segments, create one common array of segments: */ - segments=xNew(5*numsegs); + segments=xNew(5*numsegs); count=0; for(i=0;i* segment=NULL; + int numsegs; + + /*intermediary: */ + DataSet* segments_dataset=NULL; + double xnodes[3]; + double ynodes[3]; + + /*We don't know how many segments we are going to get, so have a dynamic container: */ + segments_dataset=new DataSet(); + + /*Go through elements, and call ElementSegmentsIntersection routine: */ + for(i=0;iSize(); + segments=xNew(5*numsegs); + for(i=0;i* segment=(Segment*)segments_dataset->GetObjectByOffset(i); + + /*x1,y1,x2,y2 then element_id: */ + *(segments+5*i+0)=segment->x1; + *(segments+5*i+1)=segment->y1; + *(segments+5*i+2)=segment->x2; + *(segments+5*i+3)=segment->y2; + *(segments+5*i+4)=(double)segment->eid; + } + + /*Free ressources:*/ + delete segments_dataset; + + /*Assign output pointers:*/ + *psegments=segments; + *pnumsegs=numsegs; +}/*}}}*/ + +/*Utilities*/ +void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes){/*{{{*/ + + double xsegment[2]; + double ysegment[2]; + + /*Loop through contour: */ + for(int i=0;iAddObject(new Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1])); + } + else{ + /*the segment intersected at the vertex, do not bother with this "0" length segment!:*/ + } + } + else if( (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum) ){ + + /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element, + * this will decide the coordinate: */ + if (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])){ + coord1=0; + if(edge1==IntersectEnum){coord2=alpha1;} + if(edge2==IntersectEnum){coord2=beta1;} + if(edge3==IntersectEnum){coord2=gamma1;} + } + else{ + if(edge1==IntersectEnum){coord1=alpha1;} + if(edge2==IntersectEnum){coord1=beta1;} + if(edge3==IntersectEnum){coord1=gamma1;} + coord2=1.0; + } + + xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]); + xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]); + yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]); + yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]); + + segments_dataset->AddObject(new Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1])); + } + else{ + /*No interesections, but the segment might be entirely inside this triangle!: */ + if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){ + segments_dataset->AddObject(new Segment(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1])); + } + } +}/*}}}*/ +bool NodeInElement(double* xnodes, double* ynodes, double x, double y){/*{{{*/ + + double x1,y1; + double x2,y2; + double x3,y3; + double lambda1,lambda2,lambda3; + double det; + + x1=xnodes[0]; + x2=xnodes[1]; + x3=xnodes[2]; + y1=ynodes[0]; + y2=ynodes[1]; + y3=ynodes[2]; + + /*compute determinant: */ + det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1; + + /*area coordinates: */ + lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det; + lambda2=((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/det; + lambda3=1-lambda1-lambda2; + + if( ((lambda1<=1) && (lambda1>=0)) && ((lambda2<=1) && (lambda2>=0)) && ((lambda3<=1) && (lambda3>=0)) )return true; + else return false; + +}/*}}}*/ +int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){/*{{{*/ + + /*See ISSM_DIR/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */ + + /*output: */ + double alpha=-1; + double beta=-1; + + double xA,xB,xC,xD,yA,yB,yC,yD; + double O2A[2],O2B[2],O1C[2],O1D[2]; + double n1[2],n2[2]; + double test1, test2, test3, test4; + double det; + double O2O1[2]; + double pO1A,pO1B,pO1C,pO1D; + + xA=x1[0]; yA=y1[0]; + xB=x1[1]; yB=y1[1]; + xC=x2[0]; yC=y2[0]; + xD=x2[1]; yD=y2[1]; + + O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2); + O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2); + O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2); + O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2); + + n1[0]=yA-yB; n1[1]=xB-xA; //normal vector to segA + n2[0]=yC-yD; n2[1]=xD-xC; //normal vector to segB + + test1=n2[0]*O2A[0]+n2[1]*O2A[1]; + test2=n2[0]*O2B[0]+n2[1]*O2B[1]; + + if (test1*test2>0){ + return SeparateEnum; + } + + test3=n1[0]*O1C[0]+n1[1]*O1C[1]; + test4=n1[0]*O1D[0]+n1[1]*O1D[1]; + + if (test3*test4>0){ + return SeparateEnum; + } + + /*If colinear: */ + det=n1[0]*n2[1]-n2[0]*n1[1]; + + if(test1*test2==0 && test3*test4==0 && det==0){ + + //projection on the axis O1O2 + O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2); + O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2); + + pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]); + pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]); + pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1]; + pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1]; + + //test if one point is included in the other segment (->intersects=true) + if ((pO1C-pO1A)*(pO1D-pO1A)<0){ + alpha=0; beta=0; + *palpha=alpha;*pbeta=beta; + return ColinearEnum; + } + if ((pO1C-pO1B)*(pO1D-pO1B)<0){ + alpha=0; beta=0; + *palpha=alpha;*pbeta=beta; + return ColinearEnum; + } + if ((pO1A-pO1C)*(pO1B-pO1C)<0){ + alpha=0; beta=0; + *palpha=alpha;*pbeta=beta; + return ColinearEnum; + } + if ((pO1A-pO1D)*(pO1B-pO1D)<0){ + alpha=0; beta=0; + *palpha=alpha;*pbeta=beta; + return ColinearEnum; + } + + //test if the 2 segments have the same middle (->intersects=true) + if (O2O1==0){ + alpha=0; beta=0; + *palpha=alpha;*pbeta=beta; + return ColinearEnum; + } + + //if we are here, both segments are colinear, but do not interset: + alpha=-1; beta=-1; + *palpha=alpha;*pbeta=beta; + return SeparateEnum; + } + + /*if we are here, both segments intersect. Determine where in the segment coordinate + * system: */ + beta=-1; + alpha=-(xA*yB-xC*yB+yC*xB-yC*xA+xC*yA-yA*xB)/(-xD*yB+xD*yA+xC*yB-xC*yA-yD*xA+yD*xB+yC*xA-yC*xB); //from intersect.m in formal calculus + + *palpha=alpha;*pbeta=beta; + return IntersectEnum; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h (revision 14217) +++ /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h (revision 14218) @@ -10,6 +10,6 @@ /* local prototypes: */ -void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods, Contour** contours,int numcontours); -void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes); +void MeshProfileIntersectionx(int** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods, Contour** contours,int numcontours); +void MeshSegmentsIntersection(int** 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); Index: sm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp (revision 14217) +++ (revision ) @@ -1,52 +1,0 @@ -/*! \file MeshSegmentsIntersectionx.c - */ - -#include "./MeshProfileIntersectionx.h" - -void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes){ - - int i,j; - - /*output: */ - double* segments=NULL; - Segment* segment=NULL; - int numsegs; - - /*intermediary: */ - DataSet* segments_dataset=NULL; - double xnodes[3]; - double ynodes[3]; - - /*We don't know how many segments we are going to get, so have a dynamic container: */ - segments_dataset=new DataSet(); - - /*Go through elements, and call ElementSegmentsIntersection routine: */ - for(i=0;iSize(); - segments=xNew(5*numsegs); - for(i=0;i* segment=(Segment*)segments_dataset->GetObjectByOffset(i); - - /*x1,y1,x2,y2 then element_id: */ - *(segments+5*i+0)=segment->x1; - *(segments+5*i+1)=segment->y1; - *(segments+5*i+2)=segment->x2; - *(segments+5*i+3)=segment->y2; - *(segments+5*i+4)=(double)segment->eid; - } - - /*Free ressources:*/ - delete segments_dataset; - - /*Assign output pointers:*/ - *psegments=segments; - *pnumsegs=numsegs; -} Index: sm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp (revision 14217) +++ (revision ) @@ -1,32 +1,0 @@ -/*! \file NodeInElement.cpp -*/ - -#include "./MeshProfileIntersectionx.h" - -bool NodeInElement(double* xnodes, double* ynodes, double x, double y){ - - double x1,y1; - double x2,y2; - double x3,y3; - double lambda1,lambda2,lambda3; - double det; - - x1=xnodes[0]; - x2=xnodes[1]; - x3=xnodes[2]; - y1=ynodes[0]; - y2=ynodes[1]; - y3=ynodes[2]; - - /*compute determinant: */ - det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1; - - /*area coordinates: */ - lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det; - lambda2=((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/det; - lambda3=1-lambda1-lambda2; - - if( ((lambda1<=1) && (lambda1>=0)) && ((lambda2<=1) && (lambda2>=0)) && ((lambda3<=1) && (lambda3>=0)) )return true; - else return false; - -} Index: sm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp (revision 14217) +++ (revision ) @@ -1,105 +1,0 @@ -/*! \file SegmentIntersect.cpp -*/ - -#include "./MeshProfileIntersectionx.h" - -int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){ - - /*See ISSM_DIR/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */ - - /*output: */ - double alpha=-1; - double beta=-1; - - double xA,xB,xC,xD,yA,yB,yC,yD; - double O2A[2],O2B[2],O1C[2],O1D[2]; - double n1[2],n2[2]; - double test1, test2, test3, test4; - double det; - double O2O1[2]; - double pO1A,pO1B,pO1C,pO1D; - - xA=x1[0]; yA=y1[0]; - xB=x1[1]; yB=y1[1]; - xC=x2[0]; yC=y2[0]; - xD=x2[1]; yD=y2[1]; - - O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2); - O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2); - O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2); - O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2); - - n1[0]=yA-yB; n1[1]=xB-xA; //normal vector to segA - n2[0]=yC-yD; n2[1]=xD-xC; //normal vector to segB - - test1=n2[0]*O2A[0]+n2[1]*O2A[1]; - test2=n2[0]*O2B[0]+n2[1]*O2B[1]; - - if (test1*test2>0){ - return SeparateEnum; - } - - test3=n1[0]*O1C[0]+n1[1]*O1C[1]; - test4=n1[0]*O1D[0]+n1[1]*O1D[1]; - - if (test3*test4>0){ - return SeparateEnum; - } - - /*If colinear: */ - det=n1[0]*n2[1]-n2[0]*n1[1]; - - if(test1*test2==0 && test3*test4==0 && det==0){ - - //projection on the axis O1O2 - O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2); - O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2); - - pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]); - pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]); - pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1]; - pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1]; - - //test if one point is included in the other segment (->intersects=true) - if ((pO1C-pO1A)*(pO1D-pO1A)<0){ - alpha=0; beta=0; - *palpha=alpha;*pbeta=beta; - return ColinearEnum; - } - if ((pO1C-pO1B)*(pO1D-pO1B)<0){ - alpha=0; beta=0; - *palpha=alpha;*pbeta=beta; - return ColinearEnum; - } - if ((pO1A-pO1C)*(pO1B-pO1C)<0){ - alpha=0; beta=0; - *palpha=alpha;*pbeta=beta; - return ColinearEnum; - } - if ((pO1A-pO1D)*(pO1B-pO1D)<0){ - alpha=0; beta=0; - *palpha=alpha;*pbeta=beta; - return ColinearEnum; - } - - //test if the 2 segments have the same middle (->intersects=true) - if (O2O1==0){ - alpha=0; beta=0; - *palpha=alpha;*pbeta=beta; - return ColinearEnum; - } - - //if we are here, both segments are colinear, but do not interset: - alpha=-1; beta=-1; - *palpha=alpha;*pbeta=beta; - return SeparateEnum; - } - - /*if we are here, both segments intersect. Determine where in the segment coordinate - * system: */ - beta=-1; - alpha=-(xA*yB-xC*yB+yC*xB-yC*xA+xC*yA-yA*xB)/(-xD*yB+xD*yA+xC*yB-xC*yA-yD*xA+yD*xB+yC*xA-yC*xB); //from intersect.m in formal calculus - - *palpha=alpha;*pbeta=beta; - return IntersectEnum; -} Index: sm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/intersect.m =================================================================== --- /issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/intersect.m (revision 14217) +++ (revision ) @@ -1,12 +1,0 @@ -syms xA yA xB yB xC yC xD yD alpha beta x y - -A=[xA;yA]; -B=[xB;yB]; -C=[xC;yC]; -D=[xD;yD]; - - -Eq=C+alpha*(D-C)-A+beta*(B-A); - -%from Eq, we specify the system to solve: -S=solve( xC+alpha*(xD-xC)-xA+beta*(xB-xA), yC+alpha*(yD-yC)-yA+beta*(yB-yA),alpha,beta);