Changeset 14218


Ignore:
Timestamp:
01/09/13 08:37:02 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: MeshProfileIntersection now writes integers + lots of clean up

Location:
issm/trunk-jpl/src/c
Files:
6 deleted
3 edited

Legend:

Unmodified
Added
Removed

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

    r14191 r14218  
    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\

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

    r13762 r14218  
    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;
    … …  
    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
    … …  
    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++){
    … …  
    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}/*}}}*/

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

    r13623 r14218  
    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);
Note: See TracChangeset for help on using the changeset viewer.