Changeset 4785


Ignore:
Timestamp:
07/23/10 18:41:49 (15 years ago)
Author:
Eric.Larour
Message:

Almost finished MeshProfileIntersectionx module

Location:
issm/trunk/src
Files:
6 added
6 edited

Legend:

Unmodified
Added
Removed

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

    r4773 r4785  
    506506                                        ./modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp\
    507507                                        ./modules/MeshProfileIntersectionx/ElementSegment.cpp\
     508                                        ./modules/MeshProfileIntersectionx/SegmentIntersect.cpp\
     509                                        ./modules/MeshProfileIntersectionx/NodeInElement.cpp\
    508510                                        ./modules/ContourToMeshx/ContourToMeshx.cpp\
    509511                                        ./modules/ContourToMeshx/ContourToMeshxt.cpp\
    … …  
    10651067                                        ./modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp\
    10661068                                        ./modules/MeshProfileIntersectionx/ElementSegment.cpp\
     1069                                        ./modules/MeshProfileIntersectionx/SegmentIntersect.cpp\
     1070                                        ./modules/MeshProfileIntersectionx/NodeInElement.cpp\
    10671071                                        ./modules/ContourToMeshx/ContourToMeshx.cpp\
    10681072                                        ./modules/ContourToMeshx/ContourToMeshxt.cpp\

  • issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp

    r4773 r4785  
    44#include "./MeshProfileIntersectionx.h"
    55               
    6 void ElementSegment(DataSet* segments_dataset,double* xgrids,double* ygrids,double* xsegment,double* ysegment){
    7 
     6void ElementSegment(DataSet* segments_dataset,int el,double* xgrids,double* ygrids,double* xsegment,double* ysegment){
    87
    98        /*We have a tria element (xgrids,ygrids) and a segment (xsegment,ysegment). Find whether they intersect.
    109         * If they do, create a Segment object with the intersection, and add to segments_dataset dataset: */
    1110
    12         ISSMERROR("not supported yet!");
     11        int i;
     12        double alpha;
     13        double alpha1,alpha2;
     14        double beta1,beta2;
     15        double gamma1,gamma2;
     16       
     17        int    edge1,edge2,edge3;
    1318
     19        double xel[2],yel[2];
     20        double coord1,coord2;
     21        double xfinal[2],yfinal[2];
     22
     23       
     24        /*edge 1: */
     25        xel[0]=xgrids[0];  yel[0]=ygrids[0]; xel[1]=xgrids[1];  yel[1]=ygrids[1];
     26        edge1=SegmentIntersect(&alpha1,&alpha2, xel,yel,xsegment,ysegment);
     27
     28        /*edge 2: */
     29        xel[0]=xgrids[1];  yel[0]=ygrids[1]; xel[1]=xgrids[2];  yel[1]=ygrids[2];
     30        edge2=SegmentIntersect(&beta1,&beta2, xel,yel,xsegment,ysegment);
     31
     32        /*edge 3: */
     33        xel[0]=xgrids[2];  yel[0]=ygrids[2]; xel[1]=xgrids[0];  yel[1]=ygrids[0];
     34        edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
     35
     36        if(    (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum)   ){
     37                /*This case is impossible: */
     38                ISSMERROR(" error: a line cannot go through 3 different vertices!");
     39        }
     40        else if(    ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum))   ){
     41       
     42                /*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */
     43                if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);}
     44                if((edge2==IntersectEnum) && (edge3==IntersectEnum)) {coord1=min(beta1,gamma1); coord2=max(beta1,gamma1);}
     45                if((edge3==IntersectEnum) && (edge1==IntersectEnum)) {coord1=min(gamma1,alpha1); coord2=max(gamma1,alpha1);}
     46
     47                /*check this segment did not intersect at a vertex of the tria: */
     48                if(coord1!=coord2){
     49
     50                        xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
     51                        xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
     52                        yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
     53                        yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
     54
     55                        segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
     56                }
     57        }
     58        else if(  (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)   ){
     59
     60                /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element,
     61                 * this will decide the coordinate: */
     62                if (NodeInElement(xgrids,ygrids,xsegment[0],ysegment[0])){
     63                        coord1=0;
     64                        if(edge1==IntersectEnum){coord2=alpha1;}
     65                        if(edge2==IntersectEnum){coord2=beta1;}
     66                        if(edge3==IntersectEnum){coord2=gamma1;}
     67                }
     68                else{
     69                        if(edge1==IntersectEnum){coord1=alpha1;}
     70                        if(edge2==IntersectEnum){coord1=beta1;}
     71                        if(edge3==IntersectEnum){coord1=gamma1;}
     72                        coord2=1.0;
     73                }
     74               
     75                xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
     76                xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
     77                yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
     78                yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
     79
     80                segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
     81        }
    1482}

  • issm/trunk/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp

    r4773 r4785  
    44#include "./MeshProfileIntersectionx.h"
    55               
    6 void ElementSegmentsIntersection(DataSet* segments_dataset,double* xgrids,double* ygrids,double* xc,double* yc,int numgrids){
     6void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xgrids,double* ygrids,double* xc,double* yc,int numgrids){
    77
    88        int i;
    … …  
    1818                ysegment[1]=yc[i+1];
    1919
    20                 ElementSegment(segments_dataset,xgrids,ygrids,xsegment,ysegment);
     20                ElementSegment(segments_dataset,el, xgrids,ygrids,xsegment,ysegment);
    2121
    2222        }

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

    r4773 r4785  
    1313void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour** contours,int numcontours);
    1414void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numgrids);
    15 void ElementSegmentsIntersection(DataSet* segments_dataset,double* xgrids,double* ygrids,double* xc,double* yc,int numgrids);
    16 void ElementSegment(DataSet* segments_dataset,double* xgrids,double* ygrids,double* xsegment,double* ysegment);
     15void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xgrids,double* ygrids,double* xc,double* yc,int numgrids);
     16void ElementSegment(DataSet* segments_dataset,int el,double* xgrids,double* ygrids,double* xsegment,double* ysegment);
     17int  SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2);
     18bool NodeInElement(double* xgrids, double* ygrids, double x, double y);
    1719
    1820#endif /* _MESHPROFILEINTERSECTIONX_H */

  • issm/trunk/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp

    r4773 r4785  
    2727                        ygrids[j]=y[*(index+3*i+j)];
    2828                }
    29                 ElementSegmentsIntersection(segments_dataset,xgrids,ygrids,xc,yc,numgrids);
     29                ElementSegmentsIntersection(segments_dataset,i,xgrids,ygrids,xc,yc,numgrids);
    3030        }
    3131

  • issm/trunk/src/m/utils/Geometry/SegIntersect.m

    r1 r4785  
    44%   return 1 if the two segments intersect
    55%   seg1=[x1 y1; x2 y2]
     6%   seg2=[x1 y1; x2 y2]
    67%
    78%   Usage:
Note: See TracChangeset for help on using the changeset viewer.