source: issm/oecreview/Archive/14064-14311/ISSM-14217-14218.diff@ 14312

Last change on this file since 14312 was 14312, checked in by Mathieu Morlighem, 12 years ago

Added 14064-14311
File size: 24.7 KB

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp

     
    1 /*! \file  SegmentIntersect.cpp
    2 */
    3 
    4 #include "./MeshProfileIntersectionx.h"
    5 
    6 int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){
    7 
    8         /*See ISSM_DIR/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */
    9 
    10         /*output: */
    11         double alpha=-1;
    12         double beta=-1;
    13 
    14         double xA,xB,xC,xD,yA,yB,yC,yD;
    15         double O2A[2],O2B[2],O1C[2],O1D[2];
    16         double n1[2],n2[2];
    17         double test1, test2, test3, test4;
    18         double det;
    19         double O2O1[2];
    20         double pO1A,pO1B,pO1C,pO1D;
    21 
    22         xA=x1[0]; yA=y1[0];
    23         xB=x1[1]; yB=y1[1];
    24         xC=x2[0]; yC=y2[0];
    25         xD=x2[1]; yD=y2[1];
    26 
    27         O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2);
    28         O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2);
    29         O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2);
    30         O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2);
    31 
    32         n1[0]=yA-yB; n1[1]=xB-xA;  //normal vector to segA
    33         n2[0]=yC-yD; n2[1]=xD-xC;  //normal vector to segB
    34 
    35         test1=n2[0]*O2A[0]+n2[1]*O2A[1];
    36         test2=n2[0]*O2B[0]+n2[1]*O2B[1];
    37 
    38         if (test1*test2>0){
    39                 return SeparateEnum;
    40         }
    41 
    42         test3=n1[0]*O1C[0]+n1[1]*O1C[1];
    43         test4=n1[0]*O1D[0]+n1[1]*O1D[1];
    44 
    45         if (test3*test4>0){
    46                 return SeparateEnum;
    47         }
    48 
    49         /*If colinear: */
    50         det=n1[0]*n2[1]-n2[0]*n1[1];
    51 
    52         if(test1*test2==0 && test3*test4==0 && det==0){
    53 
    54                 //projection on the axis O1O2
    55                 O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2);
    56                 O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2);
    57 
    58                 pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]);
    59                 pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]);
    60                 pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1];
    61                 pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1];
    62 
    63                 //test if one point is included in the other segment (->intersects=true)
    64                 if ((pO1C-pO1A)*(pO1D-pO1A)<0){
    65                         alpha=0; beta=0;
    66                         *palpha=alpha;*pbeta=beta;
    67                         return ColinearEnum;
    68                 }
    69                 if ((pO1C-pO1B)*(pO1D-pO1B)<0){
    70                         alpha=0; beta=0;
    71                         *palpha=alpha;*pbeta=beta;
    72                         return ColinearEnum;
    73                 }
    74                 if ((pO1A-pO1C)*(pO1B-pO1C)<0){
    75                         alpha=0; beta=0;
    76                         *palpha=alpha;*pbeta=beta;
    77                         return ColinearEnum;
    78                 }
    79                 if ((pO1A-pO1D)*(pO1B-pO1D)<0){
    80                         alpha=0; beta=0;
    81                         *palpha=alpha;*pbeta=beta;
    82                         return ColinearEnum;
    83                 }
    84 
    85                 //test if the 2 segments have the same middle (->intersects=true)
    86                 if (O2O1==0){
    87                         alpha=0; beta=0;
    88                         *palpha=alpha;*pbeta=beta;
    89                         return ColinearEnum;
    90                 }
    91 
    92                 //if we are here, both segments are colinear, but do not interset:
    93                 alpha=-1; beta=-1;
    94                 *palpha=alpha;*pbeta=beta;
    95                 return SeparateEnum;
    96         }
    97 
    98         /*if we are here, both segments intersect. Determine where in the segment coordinate
    99          * system: */
    100         beta=-1;
    101         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
    102 
    103         *palpha=alpha;*pbeta=beta;
    104         return IntersectEnum;
    105 }

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp

     
    1 /*! \file  ElementSegmentsIntersection.cpp
    2  */
    3 
    4 #include "./MeshProfileIntersectionx.h"
    5 
    6 void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes){
    7 
    8         int i;
    9         double xsegment[2];
    10         double ysegment[2];
    11 
    12         /*Loop through contour: */
    13         for(i=0;i<numnodes-1;i++){
    14 
    15                 xsegment[0]=xc[i];
    16                 xsegment[1]=xc[i+1];
    17                 ysegment[0]=yc[i];
    18                 ysegment[1]=yc[i+1];
    19 
    20                 ElementSegment(segments_dataset,el, xnodes,ynodes,xsegment,ysegment);
    21 
    22         }
    23 }

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp

     
    1 /*! \file  NodeInElement.cpp
    2 */
    3 
    4 #include "./MeshProfileIntersectionx.h"
    5 
    6 bool NodeInElement(double* xnodes, double* ynodes, double x, double y){
    7 
    8         double x1,y1;
    9         double x2,y2;
    10         double x3,y3;
    11         double lambda1,lambda2,lambda3;
    12         double det;
    13 
    14         x1=xnodes[0];
    15         x2=xnodes[1];
    16         x3=xnodes[2];
    17         y1=ynodes[0];
    18         y2=ynodes[1];
    19         y3=ynodes[2];
    20 
    21         /*compute determinant: */
    22         det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1;
    23 
    24         /*area coordinates: */
    25         lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det;
    26         lambda2=((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/det;
    27         lambda3=1-lambda1-lambda2;
    28 
    29         if( ((lambda1<=1) && (lambda1>=0)) && ((lambda2<=1) && (lambda2>=0)) && ((lambda3<=1) && (lambda3>=0))  )return true;
    30         else return false;
    31 
    32 }

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp

     
    1 /*! \file  MeshSegmentsIntersectionx.c
    2  */
    3 
    4 #include "./MeshProfileIntersectionx.h"
    5 
    6 void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes){
    7 
    8         int      i,j;
    9 
    10         /*output: */
    11         double*  segments=NULL;
    12         Segment<double>* segment=NULL;
    13         int     numsegs;
    14 
    15         /*intermediary: */
    16         DataSet* segments_dataset=NULL;
    17         double   xnodes[3];
    18         double   ynodes[3];
    19 
    20         /*We don't know how many segments  we are going to get, so have a dynamic container: */
    21         segments_dataset=new DataSet();
    22 
    23         /*Go through elements, and call ElementSegmentsIntersection routine: */
    24         for(i=0;i<nel;i++){
    25                 for(j=0;j<3;j++){
    26                         xnodes[j]=x[*(index+3*i+j)];
    27                         ynodes[j]=y[*(index+3*i+j)];
    28                 }
    29                 ElementSegmentsIntersection(segments_dataset,i,xnodes,ynodes,xc,yc,numnodes);
    30         }
    31 
    32         /*Using the segments_dataset dataset, create segments: */
    33         numsegs=segments_dataset->Size();
    34         segments=xNew<double>(5*numsegs);
    35         for(i=0;i<numsegs;i++){
    36                 Segment<double>* segment=(Segment<double>*)segments_dataset->GetObjectByOffset(i);
    37 
    38                 /*x1,y1,x2,y2 then element_id: */
    39                 *(segments+5*i+0)=segment->x1;
    40                 *(segments+5*i+1)=segment->y1;
    41                 *(segments+5*i+2)=segment->x2;
    42                 *(segments+5*i+3)=segment->y2;
    43                 *(segments+5*i+4)=(double)segment->eid;
    44         }
    45 
    46         /*Free ressources:*/
    47         delete segments_dataset;
    48 
    49         /*Assign output pointers:*/
    50         *psegments=segments;
    51         *pnumsegs=numsegs;
    52 }

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/intersect.m

     
    1 syms xA yA xB yB xC yC xD yD alpha beta x y
    2 
    3 A=[xA;yA];
    4 B=[xB;yB];
    5 C=[xC;yC];
    6 D=[xD;yD];
    7 
    8 
    9 Eq=C+alpha*(D-C)-A+beta*(B-A);
    10 
    11 %from Eq, we specify the system to solve:
    12 S=solve( xC+alpha*(xD-xC)-xA+beta*(xB-xA), yC+alpha*(yD-yC)-yA+beta*(yB-yA),alpha,beta);

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp

     
    1 /*! \file  ElementSegment.cpp
    2  */
    3 
    4 #include "./MeshProfileIntersectionx.h"
    5 
    6 void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment){
    7 
    8         /*We have a tria element (xnodes,ynodes) and a segment (xsegment,ysegment). Find whether they intersect.
    9          * If they do, create a Segment object with the intersection, and add to segments_dataset dataset: */
    10 
    11         double alpha1,alpha2;
    12         double beta1,beta2;
    13         double gamma1,gamma2;
    14 
    15         int    edge1,edge2,edge3;
    16 
    17         double xel[2],yel[2];
    18         double coord1,coord2;
    19         double xfinal[2],yfinal[2];
    20 
    21         /*edge 1: */
    22         xel[0]=xnodes[0];  yel[0]=ynodes[0]; xel[1]=xnodes[1];  yel[1]=ynodes[1];
    23         edge1=SegmentIntersect(&alpha1,&alpha2, xel,yel,xsegment,ysegment); //alpha1: segment coordinate of intersection. alpha2: same thing for second interesection if it exists (colinear edges)
    24 
    25         /*edge 2: */
    26         xel[0]=xnodes[1];  yel[0]=ynodes[1]; xel[1]=xnodes[2];  yel[1]=ynodes[2];
    27         edge2=SegmentIntersect(&beta1,&beta2, xel,yel,xsegment,ysegment);
    28 
    29         /*edge 3: */
    30         xel[0]=xnodes[2];  yel[0]=ynodes[2]; xel[1]=xnodes[0];  yel[1]=ynodes[0];
    31         edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
    32 
    33         /*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): */
    34 
    35         if(    (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum)   ){
    36                 /*This case is impossible: */
    37                 _error_("error: a line cannot go through 3 different vertices!");
    38         }
    39         else if(    ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum))   ){
    40 
    41                 /*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */
    42                 if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);}
    43                 if((edge2==IntersectEnum) && (edge3==IntersectEnum)) {coord1=min(beta1,gamma1); coord2=max(beta1,gamma1);}
    44                 if((edge3==IntersectEnum) && (edge1==IntersectEnum)) {coord1=min(gamma1,alpha1); coord2=max(gamma1,alpha1);}
    45 
    46                 /*check this segment did not intersect at a vertex of the tria: */
    47                 if(coord1!=coord2){
    48 
    49                         xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
    50                         xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
    51                         yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
    52                         yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
    53 
    54                         segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
    55                 }
    56                 else{
    57                         /*the segment intersected at the vertex, do not bother with this "0" length segment!:*/
    58                 }
    59         }
    60         else if(  (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)   ){
    61 
    62                 /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element,
    63                  * this will decide the coordinate: */
    64                 if (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])){
    65                         coord1=0;
    66                         if(edge1==IntersectEnum){coord2=alpha1;}
    67                         if(edge2==IntersectEnum){coord2=beta1;}
    68                         if(edge3==IntersectEnum){coord2=gamma1;}
    69                 }
    70                 else{
    71                         if(edge1==IntersectEnum){coord1=alpha1;}
    72                         if(edge2==IntersectEnum){coord1=beta1;}
    73                         if(edge3==IntersectEnum){coord1=gamma1;}
    74                         coord2=1.0;
    75                 }
    76 
    77                 xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
    78                 xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
    79                 yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
    80                 yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
    81 
    82                 segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
    83         }
    84         else{
    85                 /*No interesections, but the segment might be entirely inside this triangle!: */
    86                 if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){
    87                         segments_dataset->AddObject(new  Segment<double>(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
    88                 }
    89         }
    90 }

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h

     
    99#include "../../classes/objects/objects.h"
    1010
    1111/* local prototypes: */
    12 void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours);
    13 void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes);
     12void MeshProfileIntersectionx(int** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours);
     13void MeshSegmentsIntersection(int** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes);
    1414void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes);
    1515void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment);
    1616int  SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2);

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp

     
    33
    44#include "./MeshProfileIntersectionx.h"
    55
    6 void MeshProfileIntersectionx( double** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours){
     6void MeshProfileIntersectionx(int** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours){/*{{{*/
    77
    88        int i,j,k;
    99
    … …  
    1414        double*  yc=NULL;
    1515
    1616        /*output: */
    17         double* segments=NULL;
    18         int     numsegs;
     17        int* segments=NULL;
     18        int  numsegs;
    1919
    2020        /*intermediary: */
    21         double** allsegments=NULL;
    22         double* segmentsi=NULL;
    23         int*    allnumsegs=NULL;
    24         int     numsegsi;
    25         int     count;
     21        int** allsegments=NULL;
     22        int* segmentsi=NULL;
     23        int*  allnumsegs=NULL;
     24        int   numsegsi;
     25        int   count;
    2626
    2727        /*Allocate: */
    28         allsegments=xNew<double*>(numcontours);
     28        allsegments=xNew<int*>(numcontours);
    2929        allnumsegs=xNew<int>(numcontours);
    3030
    3131        /*Loop through all contours: */
    … …  
    5050        for(i=0;i<numcontours;i++)numsegs+=allnumsegs[i];
    5151
    5252        /*Out of all segments, create one common array of segments: */
    53         segments=xNew<double>(5*numsegs);
     53        segments=xNew<int>(5*numsegs);
    5454        count=0;
    5555        for(i=0;i<numcontours;i++){
    5656
    … …  
    6868        /*Assign output pointers:*/
    6969        *psegments=segments;
    7070        *pnumsegs=numsegs;
    71 }
     71}/*}}}*/
     72void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes){/*{{{*/
     73
     74        int      i,j;
     75
     76        /*output: */
     77        double*  segments=NULL;
     78        Segment<double>* segment=NULL;
     79        int     numsegs;
     80
     81        /*intermediary: */
     82        DataSet* segments_dataset=NULL;
     83        double   xnodes[3];
     84        double   ynodes[3];
     85
     86        /*We don't know how many segments  we are going to get, so have a dynamic container: */
     87        segments_dataset=new DataSet();
     88
     89        /*Go through elements, and call ElementSegmentsIntersection routine: */
     90        for(i=0;i<nel;i++){
     91                for(j=0;j<3;j++){
     92                        xnodes[j]=x[*(index+3*i+j)];
     93                        ynodes[j]=y[*(index+3*i+j)];
     94                }
     95                ElementSegmentsIntersection(segments_dataset,i,xnodes,ynodes,xc,yc,numnodes);
     96        }
     97
     98        /*Using the segments_dataset dataset, create segments: */
     99        numsegs=segments_dataset->Size();
     100        segments=xNew<double>(5*numsegs);
     101        for(i=0;i<numsegs;i++){
     102                Segment<double>* segment=(Segment<double>*)segments_dataset->GetObjectByOffset(i);
     103
     104                /*x1,y1,x2,y2 then element_id: */
     105                *(segments+5*i+0)=segment->x1;
     106                *(segments+5*i+1)=segment->y1;
     107                *(segments+5*i+2)=segment->x2;
     108                *(segments+5*i+3)=segment->y2;
     109                *(segments+5*i+4)=(double)segment->eid;
     110        }
     111
     112        /*Free ressources:*/
     113        delete segments_dataset;
     114
     115        /*Assign output pointers:*/
     116        *psegments=segments;
     117        *pnumsegs=numsegs;
     118}/*}}}*/
     119
     120/*Utilities*/
     121void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes){/*{{{*/
     122
     123        double xsegment[2];
     124        double ysegment[2];
     125
     126        /*Loop through contour: */
     127        for(int i=0;i<numnodes-1;i++){
     128                xsegment[0]=xc[i];
     129                xsegment[1]=xc[i+1];
     130                ysegment[0]=yc[i];
     131                ysegment[1]=yc[i+1];
     132                ElementSegment(segments_dataset,el, xnodes,ynodes,xsegment,ysegment);
     133        }
     134}/*}}}*/
     135void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment){/*{{{*/
     136
     137        /*We have a tria element (xnodes,ynodes) and a segment (xsegment,ysegment). Find whether they intersect.
     138         * If they do, create a Segment object with the intersection, and add to segments_dataset dataset: */
     139
     140        double alpha1,alpha2;
     141        double beta1,beta2;
     142        double gamma1,gamma2;
     143
     144        int    edge1,edge2,edge3;
     145
     146        double xel[2],yel[2];
     147        double coord1,coord2;
     148        double xfinal[2],yfinal[2];
     149
     150        /*edge 1: */
     151        xel[0]=xnodes[0];  yel[0]=ynodes[0]; xel[1]=xnodes[1];  yel[1]=ynodes[1];
     152        edge1=SegmentIntersect(&alpha1,&alpha2, xel,yel,xsegment,ysegment); //alpha1: segment coordinate of intersection. alpha2: same thing for second interesection if it exists (colinear edges)
     153
     154        /*edge 2: */
     155        xel[0]=xnodes[1];  yel[0]=ynodes[1]; xel[1]=xnodes[2];  yel[1]=ynodes[2];
     156        edge2=SegmentIntersect(&beta1,&beta2, xel,yel,xsegment,ysegment);
     157
     158        /*edge 3: */
     159        xel[0]=xnodes[2];  yel[0]=ynodes[2]; xel[1]=xnodes[0];  yel[1]=ynodes[0];
     160        edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
     161
     162        /*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): */
     163
     164        if(    (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum)   ){
     165                /*This case is impossible: */
     166                _error_("error: a line cannot go through 3 different vertices!");
     167        }
     168        else if(    ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum))   ){
     169
     170                /*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */
     171                if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);}
     172                if((edge2==IntersectEnum) && (edge3==IntersectEnum)) {coord1=min(beta1,gamma1); coord2=max(beta1,gamma1);}
     173                if((edge3==IntersectEnum) && (edge1==IntersectEnum)) {coord1=min(gamma1,alpha1); coord2=max(gamma1,alpha1);}
     174
     175                /*check this segment did not intersect at a vertex of the tria: */
     176                if(coord1!=coord2){
     177
     178                        xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
     179                        xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
     180                        yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
     181                        yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
     182
     183                        segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
     184                }
     185                else{
     186                        /*the segment intersected at the vertex, do not bother with this "0" length segment!:*/
     187                }
     188        }
     189        else if(  (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)   ){
     190
     191                /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element,
     192                 * this will decide the coordinate: */
     193                if (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])){
     194                        coord1=0;
     195                        if(edge1==IntersectEnum){coord2=alpha1;}
     196                        if(edge2==IntersectEnum){coord2=beta1;}
     197                        if(edge3==IntersectEnum){coord2=gamma1;}
     198                }
     199                else{
     200                        if(edge1==IntersectEnum){coord1=alpha1;}
     201                        if(edge2==IntersectEnum){coord1=beta1;}
     202                        if(edge3==IntersectEnum){coord1=gamma1;}
     203                        coord2=1.0;
     204                }
     205
     206                xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
     207                xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
     208                yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
     209                yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
     210
     211                segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
     212        }
     213        else{
     214                /*No interesections, but the segment might be entirely inside this triangle!: */
     215                if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){
     216                        segments_dataset->AddObject(new  Segment<double>(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
     217                }
     218        }
     219}/*}}}*/
     220bool NodeInElement(double* xnodes, double* ynodes, double x, double y){/*{{{*/
     221
     222        double x1,y1;
     223        double x2,y2;
     224        double x3,y3;
     225        double lambda1,lambda2,lambda3;
     226        double det;
     227
     228        x1=xnodes[0];
     229        x2=xnodes[1];
     230        x3=xnodes[2];
     231        y1=ynodes[0];
     232        y2=ynodes[1];
     233        y3=ynodes[2];
     234
     235        /*compute determinant: */
     236        det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1;
     237
     238        /*area coordinates: */
     239        lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det;
     240        lambda2=((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/det;
     241        lambda3=1-lambda1-lambda2;
     242
     243        if( ((lambda1<=1) && (lambda1>=0)) && ((lambda2<=1) && (lambda2>=0)) && ((lambda3<=1) && (lambda3>=0))  )return true;
     244        else return false;
     245
     246}/*}}}*/
     247int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){/*{{{*/
     248
     249        /*See ISSM_DIR/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */
     250
     251        /*output: */
     252        double alpha=-1;
     253        double beta=-1;
     254
     255        double xA,xB,xC,xD,yA,yB,yC,yD;
     256        double O2A[2],O2B[2],O1C[2],O1D[2];
     257        double n1[2],n2[2];
     258        double test1, test2, test3, test4;
     259        double det;
     260        double O2O1[2];
     261        double pO1A,pO1B,pO1C,pO1D;
     262
     263        xA=x1[0]; yA=y1[0];
     264        xB=x1[1]; yB=y1[1];
     265        xC=x2[0]; yC=y2[0];
     266        xD=x2[1]; yD=y2[1];
     267
     268        O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2);
     269        O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2);
     270        O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2);
     271        O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2);
     272
     273        n1[0]=yA-yB; n1[1]=xB-xA;  //normal vector to segA
     274        n2[0]=yC-yD; n2[1]=xD-xC;  //normal vector to segB
     275
     276        test1=n2[0]*O2A[0]+n2[1]*O2A[1];
     277        test2=n2[0]*O2B[0]+n2[1]*O2B[1];
     278
     279        if (test1*test2>0){
     280                return SeparateEnum;
     281        }
     282
     283        test3=n1[0]*O1C[0]+n1[1]*O1C[1];
     284        test4=n1[0]*O1D[0]+n1[1]*O1D[1];
     285
     286        if (test3*test4>0){
     287                return SeparateEnum;
     288        }
     289
     290        /*If colinear: */
     291        det=n1[0]*n2[1]-n2[0]*n1[1];
     292
     293        if(test1*test2==0 && test3*test4==0 && det==0){
     294
     295                //projection on the axis O1O2
     296                O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2);
     297                O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2);
     298
     299                pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]);
     300                pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]);
     301                pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1];
     302                pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1];
     303
     304                //test if one point is included in the other segment (->intersects=true)
     305                if ((pO1C-pO1A)*(pO1D-pO1A)<0){
     306                        alpha=0; beta=0;
     307                        *palpha=alpha;*pbeta=beta;
     308                        return ColinearEnum;
     309                }
     310                if ((pO1C-pO1B)*(pO1D-pO1B)<0){
     311                        alpha=0; beta=0;
     312                        *palpha=alpha;*pbeta=beta;
     313                        return ColinearEnum;
     314                }
     315                if ((pO1A-pO1C)*(pO1B-pO1C)<0){
     316                        alpha=0; beta=0;
     317                        *palpha=alpha;*pbeta=beta;
     318                        return ColinearEnum;
     319                }
     320                if ((pO1A-pO1D)*(pO1B-pO1D)<0){
     321                        alpha=0; beta=0;
     322                        *palpha=alpha;*pbeta=beta;
     323                        return ColinearEnum;
     324                }
     325
     326                //test if the 2 segments have the same middle (->intersects=true)
     327                if (O2O1==0){
     328                        alpha=0; beta=0;
     329                        *palpha=alpha;*pbeta=beta;
     330                        return ColinearEnum;
     331                }
     332
     333                //if we are here, both segments are colinear, but do not interset:
     334                alpha=-1; beta=-1;
     335                *palpha=alpha;*pbeta=beta;
     336                return SeparateEnum;
     337        }
     338
     339        /*if we are here, both segments intersect. Determine where in the segment coordinate
     340         * system: */
     341        beta=-1;
     342        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
     343
     344        *palpha=alpha;*pbeta=beta;
     345        return IntersectEnum;
     346}/*}}}*/

  • ../trunk-jpl/src/c/Makefile.am

     
    780780                        ./modules/AverageFilterx/AverageFilterx.h\
    781781                        ./modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp\
    782782                        ./modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h\
    783                         ./modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp\
    784                         ./modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp\
    785                         ./modules/MeshProfileIntersectionx/ElementSegment.cpp\
    786                         ./modules/MeshProfileIntersectionx/SegmentIntersect.cpp\
    787                         ./modules/MeshProfileIntersectionx/NodeInElement.cpp\
    788783                        ./modules/ContourToMeshx/ContourToMeshx.cpp\
    789784                        ./modules/ContourToMeshx/ContourToMeshxt.cpp\
    790785                        ./modules/ContourToMeshx/ContourToMeshx.h\
Note: See TracBrowser for help on using the repository browser.