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
RevLine 
[14312]1Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp
2===================================================================
3--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp (revision 14217)
4+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/SegmentIntersect.cpp (revision 14218)
5@@ -1,105 +0,0 @@
6-/*! \file SegmentIntersect.cpp
7-*/
8-
9-#include "./MeshProfileIntersectionx.h"
10-
11-int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){
12-
13- /*See ISSM_DIR/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */
14-
15- /*output: */
16- double alpha=-1;
17- double beta=-1;
18-
19- double xA,xB,xC,xD,yA,yB,yC,yD;
20- double O2A[2],O2B[2],O1C[2],O1D[2];
21- double n1[2],n2[2];
22- double test1, test2, test3, test4;
23- double det;
24- double O2O1[2];
25- double pO1A,pO1B,pO1C,pO1D;
26-
27- xA=x1[0]; yA=y1[0];
28- xB=x1[1]; yB=y1[1];
29- xC=x2[0]; yC=y2[0];
30- xD=x2[1]; yD=y2[1];
31-
32- O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2);
33- O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2);
34- O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2);
35- O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2);
36-
37- n1[0]=yA-yB; n1[1]=xB-xA; //normal vector to segA
38- n2[0]=yC-yD; n2[1]=xD-xC; //normal vector to segB
39-
40- test1=n2[0]*O2A[0]+n2[1]*O2A[1];
41- test2=n2[0]*O2B[0]+n2[1]*O2B[1];
42-
43- if (test1*test2>0){
44- return SeparateEnum;
45- }
46-
47- test3=n1[0]*O1C[0]+n1[1]*O1C[1];
48- test4=n1[0]*O1D[0]+n1[1]*O1D[1];
49-
50- if (test3*test4>0){
51- return SeparateEnum;
52- }
53-
54- /*If colinear: */
55- det=n1[0]*n2[1]-n2[0]*n1[1];
56-
57- if(test1*test2==0 && test3*test4==0 && det==0){
58-
59- //projection on the axis O1O2
60- O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2);
61- O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2);
62-
63- pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]);
64- pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]);
65- pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1];
66- pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1];
67-
68- //test if one point is included in the other segment (->intersects=true)
69- if ((pO1C-pO1A)*(pO1D-pO1A)<0){
70- alpha=0; beta=0;
71- *palpha=alpha;*pbeta=beta;
72- return ColinearEnum;
73- }
74- if ((pO1C-pO1B)*(pO1D-pO1B)<0){
75- alpha=0; beta=0;
76- *palpha=alpha;*pbeta=beta;
77- return ColinearEnum;
78- }
79- if ((pO1A-pO1C)*(pO1B-pO1C)<0){
80- alpha=0; beta=0;
81- *palpha=alpha;*pbeta=beta;
82- return ColinearEnum;
83- }
84- if ((pO1A-pO1D)*(pO1B-pO1D)<0){
85- alpha=0; beta=0;
86- *palpha=alpha;*pbeta=beta;
87- return ColinearEnum;
88- }
89-
90- //test if the 2 segments have the same middle (->intersects=true)
91- if (O2O1==0){
92- alpha=0; beta=0;
93- *palpha=alpha;*pbeta=beta;
94- return ColinearEnum;
95- }
96-
97- //if we are here, both segments are colinear, but do not interset:
98- alpha=-1; beta=-1;
99- *palpha=alpha;*pbeta=beta;
100- return SeparateEnum;
101- }
102-
103- /*if we are here, both segments intersect. Determine where in the segment coordinate
104- * system: */
105- beta=-1;
106- 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
107-
108- *palpha=alpha;*pbeta=beta;
109- return IntersectEnum;
110-}
111Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp
112===================================================================
113--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp (revision 14217)
114+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp (revision 14218)
115@@ -1,23 +0,0 @@
116-/*! \file ElementSegmentsIntersection.cpp
117- */
118-
119-#include "./MeshProfileIntersectionx.h"
120-
121-void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes){
122-
123- int i;
124- double xsegment[2];
125- double ysegment[2];
126-
127- /*Loop through contour: */
128- for(i=0;i<numnodes-1;i++){
129-
130- xsegment[0]=xc[i];
131- xsegment[1]=xc[i+1];
132- ysegment[0]=yc[i];
133- ysegment[1]=yc[i+1];
134-
135- ElementSegment(segments_dataset,el, xnodes,ynodes,xsegment,ysegment);
136-
137- }
138-}
139Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp
140===================================================================
141--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp (revision 14217)
142+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/NodeInElement.cpp (revision 14218)
143@@ -1,32 +0,0 @@
144-/*! \file NodeInElement.cpp
145-*/
146-
147-#include "./MeshProfileIntersectionx.h"
148-
149-bool NodeInElement(double* xnodes, double* ynodes, double x, double y){
150-
151- double x1,y1;
152- double x2,y2;
153- double x3,y3;
154- double lambda1,lambda2,lambda3;
155- double det;
156-
157- x1=xnodes[0];
158- x2=xnodes[1];
159- x3=xnodes[2];
160- y1=ynodes[0];
161- y2=ynodes[1];
162- y3=ynodes[2];
163-
164- /*compute determinant: */
165- det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1;
166-
167- /*area coordinates: */
168- lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det;
169- lambda2=((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/det;
170- lambda3=1-lambda1-lambda2;
171-
172- if( ((lambda1<=1) && (lambda1>=0)) && ((lambda2<=1) && (lambda2>=0)) && ((lambda3<=1) && (lambda3>=0)) )return true;
173- else return false;
174-
175-}
176Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp
177===================================================================
178--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp (revision 14217)
179+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp (revision 14218)
180@@ -1,52 +0,0 @@
181-/*! \file MeshSegmentsIntersectionx.c
182- */
183-
184-#include "./MeshProfileIntersectionx.h"
185-
186-void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes){
187-
188- int i,j;
189-
190- /*output: */
191- double* segments=NULL;
192- Segment<double>* segment=NULL;
193- int numsegs;
194-
195- /*intermediary: */
196- DataSet* segments_dataset=NULL;
197- double xnodes[3];
198- double ynodes[3];
199-
200- /*We don't know how many segments we are going to get, so have a dynamic container: */
201- segments_dataset=new DataSet();
202-
203- /*Go through elements, and call ElementSegmentsIntersection routine: */
204- for(i=0;i<nel;i++){
205- for(j=0;j<3;j++){
206- xnodes[j]=x[*(index+3*i+j)];
207- ynodes[j]=y[*(index+3*i+j)];
208- }
209- ElementSegmentsIntersection(segments_dataset,i,xnodes,ynodes,xc,yc,numnodes);
210- }
211-
212- /*Using the segments_dataset dataset, create segments: */
213- numsegs=segments_dataset->Size();
214- segments=xNew<double>(5*numsegs);
215- for(i=0;i<numsegs;i++){
216- Segment<double>* segment=(Segment<double>*)segments_dataset->GetObjectByOffset(i);
217-
218- /*x1,y1,x2,y2 then element_id: */
219- *(segments+5*i+0)=segment->x1;
220- *(segments+5*i+1)=segment->y1;
221- *(segments+5*i+2)=segment->x2;
222- *(segments+5*i+3)=segment->y2;
223- *(segments+5*i+4)=(double)segment->eid;
224- }
225-
226- /*Free ressources:*/
227- delete segments_dataset;
228-
229- /*Assign output pointers:*/
230- *psegments=segments;
231- *pnumsegs=numsegs;
232-}
233Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/intersect.m
234===================================================================
235--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/intersect.m (revision 14217)
236+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/intersect.m (revision 14218)
237@@ -1,12 +0,0 @@
238-syms xA yA xB yB xC yC xD yD alpha beta x y
239-
240-A=[xA;yA];
241-B=[xB;yB];
242-C=[xC;yC];
243-D=[xD;yD];
244-
245-
246-Eq=C+alpha*(D-C)-A+beta*(B-A);
247-
248-%from Eq, we specify the system to solve:
249-S=solve( xC+alpha*(xD-xC)-xA+beta*(xB-xA), yC+alpha*(yD-yC)-yA+beta*(yB-yA),alpha,beta);
250Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp
251===================================================================
252--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp (revision 14217)
253+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp (revision 14218)
254@@ -1,90 +0,0 @@
255-/*! \file ElementSegment.cpp
256- */
257-
258-#include "./MeshProfileIntersectionx.h"
259-
260-void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment){
261-
262- /*We have a tria element (xnodes,ynodes) and a segment (xsegment,ysegment). Find whether they intersect.
263- * If they do, create a Segment object with the intersection, and add to segments_dataset dataset: */
264-
265- double alpha1,alpha2;
266- double beta1,beta2;
267- double gamma1,gamma2;
268-
269- int edge1,edge2,edge3;
270-
271- double xel[2],yel[2];
272- double coord1,coord2;
273- double xfinal[2],yfinal[2];
274-
275- /*edge 1: */
276- xel[0]=xnodes[0]; yel[0]=ynodes[0]; xel[1]=xnodes[1]; yel[1]=ynodes[1];
277- edge1=SegmentIntersect(&alpha1,&alpha2, xel,yel,xsegment,ysegment); //alpha1: segment coordinate of intersection. alpha2: same thing for second interesection if it exists (colinear edges)
278-
279- /*edge 2: */
280- xel[0]=xnodes[1]; yel[0]=ynodes[1]; xel[1]=xnodes[2]; yel[1]=ynodes[2];
281- edge2=SegmentIntersect(&beta1,&beta2, xel,yel,xsegment,ysegment);
282-
283- /*edge 3: */
284- xel[0]=xnodes[2]; yel[0]=ynodes[2]; xel[1]=xnodes[0]; yel[1]=ynodes[0];
285- edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
286-
287- /*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): */
288-
289- if( (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum) ){
290- /*This case is impossible: */
291- _error_("error: a line cannot go through 3 different vertices!");
292- }
293- else if( ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum)) ){
294-
295- /*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */
296- if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);}
297- if((edge2==IntersectEnum) && (edge3==IntersectEnum)) {coord1=min(beta1,gamma1); coord2=max(beta1,gamma1);}
298- if((edge3==IntersectEnum) && (edge1==IntersectEnum)) {coord1=min(gamma1,alpha1); coord2=max(gamma1,alpha1);}
299-
300- /*check this segment did not intersect at a vertex of the tria: */
301- if(coord1!=coord2){
302-
303- xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
304- xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
305- yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
306- yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
307-
308- segments_dataset->AddObject(new Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
309- }
310- else{
311- /*the segment intersected at the vertex, do not bother with this "0" length segment!:*/
312- }
313- }
314- else if( (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum) ){
315-
316- /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element,
317- * this will decide the coordinate: */
318- if (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])){
319- coord1=0;
320- if(edge1==IntersectEnum){coord2=alpha1;}
321- if(edge2==IntersectEnum){coord2=beta1;}
322- if(edge3==IntersectEnum){coord2=gamma1;}
323- }
324- else{
325- if(edge1==IntersectEnum){coord1=alpha1;}
326- if(edge2==IntersectEnum){coord1=beta1;}
327- if(edge3==IntersectEnum){coord1=gamma1;}
328- coord2=1.0;
329- }
330-
331- xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
332- xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
333- yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
334- yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
335-
336- segments_dataset->AddObject(new Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
337- }
338- else{
339- /*No interesections, but the segment might be entirely inside this triangle!: */
340- if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){
341- segments_dataset->AddObject(new Segment<double>(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
342- }
343- }
344-}
345Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h
346===================================================================
347--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h (revision 14217)
348+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h (revision 14218)
349@@ -9,8 +9,8 @@
350 #include "../../classes/objects/objects.h"
351
352 /* local prototypes: */
353-void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods, Contour<IssmPDouble>** contours,int numcontours);
354-void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes);
355+void MeshProfileIntersectionx(int** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods, Contour<IssmPDouble>** contours,int numcontours);
356+void MeshSegmentsIntersection(int** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes);
357 void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes);
358 void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment);
359 int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2);
360Index: ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp
361===================================================================
362--- ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 14217)
363+++ ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp (revision 14218)
364@@ -3,7 +3,7 @@
365
366 #include "./MeshProfileIntersectionx.h"
367
368-void MeshProfileIntersectionx( double** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods, Contour<IssmPDouble>** contours,int numcontours){
369+void MeshProfileIntersectionx(int** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods, Contour<IssmPDouble>** contours,int numcontours){/*{{{*/
370
371 int i,j,k;
372
373@@ -14,18 +14,18 @@
374 double* yc=NULL;
375
376 /*output: */
377- double* segments=NULL;
378- int numsegs;
379+ int* segments=NULL;
380+ int numsegs;
381
382 /*intermediary: */
383- double** allsegments=NULL;
384- double* segmentsi=NULL;
385- int* allnumsegs=NULL;
386- int numsegsi;
387- int count;
388+ int** allsegments=NULL;
389+ int* segmentsi=NULL;
390+ int* allnumsegs=NULL;
391+ int numsegsi;
392+ int count;
393
394 /*Allocate: */
395- allsegments=xNew<double*>(numcontours);
396+ allsegments=xNew<int*>(numcontours);
397 allnumsegs=xNew<int>(numcontours);
398
399 /*Loop through all contours: */
400@@ -50,7 +50,7 @@
401 for(i=0;i<numcontours;i++)numsegs+=allnumsegs[i];
402
403 /*Out of all segments, create one common array of segments: */
404- segments=xNew<double>(5*numsegs);
405+ segments=xNew<int>(5*numsegs);
406 count=0;
407 for(i=0;i<numcontours;i++){
408
409@@ -68,4 +68,279 @@
410 /*Assign output pointers:*/
411 *psegments=segments;
412 *pnumsegs=numsegs;
413-}
414+}/*}}}*/
415+void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes){/*{{{*/
416+
417+ int i,j;
418+
419+ /*output: */
420+ double* segments=NULL;
421+ Segment<double>* segment=NULL;
422+ int numsegs;
423+
424+ /*intermediary: */
425+ DataSet* segments_dataset=NULL;
426+ double xnodes[3];
427+ double ynodes[3];
428+
429+ /*We don't know how many segments we are going to get, so have a dynamic container: */
430+ segments_dataset=new DataSet();
431+
432+ /*Go through elements, and call ElementSegmentsIntersection routine: */
433+ for(i=0;i<nel;i++){
434+ for(j=0;j<3;j++){
435+ xnodes[j]=x[*(index+3*i+j)];
436+ ynodes[j]=y[*(index+3*i+j)];
437+ }
438+ ElementSegmentsIntersection(segments_dataset,i,xnodes,ynodes,xc,yc,numnodes);
439+ }
440+
441+ /*Using the segments_dataset dataset, create segments: */
442+ numsegs=segments_dataset->Size();
443+ segments=xNew<double>(5*numsegs);
444+ for(i=0;i<numsegs;i++){
445+ Segment<double>* segment=(Segment<double>*)segments_dataset->GetObjectByOffset(i);
446+
447+ /*x1,y1,x2,y2 then element_id: */
448+ *(segments+5*i+0)=segment->x1;
449+ *(segments+5*i+1)=segment->y1;
450+ *(segments+5*i+2)=segment->x2;
451+ *(segments+5*i+3)=segment->y2;
452+ *(segments+5*i+4)=(double)segment->eid;
453+ }
454+
455+ /*Free ressources:*/
456+ delete segments_dataset;
457+
458+ /*Assign output pointers:*/
459+ *psegments=segments;
460+ *pnumsegs=numsegs;
461+}/*}}}*/
462+
463+/*Utilities*/
464+void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes){/*{{{*/
465+
466+ double xsegment[2];
467+ double ysegment[2];
468+
469+ /*Loop through contour: */
470+ for(int i=0;i<numnodes-1;i++){
471+ xsegment[0]=xc[i];
472+ xsegment[1]=xc[i+1];
473+ ysegment[0]=yc[i];
474+ ysegment[1]=yc[i+1];
475+ ElementSegment(segments_dataset,el, xnodes,ynodes,xsegment,ysegment);
476+ }
477+}/*}}}*/
478+void ElementSegment(DataSet* segments_dataset,int el,double* xnodes,double* ynodes,double* xsegment,double* ysegment){/*{{{*/
479+
480+ /*We have a tria element (xnodes,ynodes) and a segment (xsegment,ysegment). Find whether they intersect.
481+ * If they do, create a Segment object with the intersection, and add to segments_dataset dataset: */
482+
483+ double alpha1,alpha2;
484+ double beta1,beta2;
485+ double gamma1,gamma2;
486+
487+ int edge1,edge2,edge3;
488+
489+ double xel[2],yel[2];
490+ double coord1,coord2;
491+ double xfinal[2],yfinal[2];
492+
493+ /*edge 1: */
494+ xel[0]=xnodes[0]; yel[0]=ynodes[0]; xel[1]=xnodes[1]; yel[1]=ynodes[1];
495+ edge1=SegmentIntersect(&alpha1,&alpha2, xel,yel,xsegment,ysegment); //alpha1: segment coordinate of intersection. alpha2: same thing for second interesection if it exists (colinear edges)
496+
497+ /*edge 2: */
498+ xel[0]=xnodes[1]; yel[0]=ynodes[1]; xel[1]=xnodes[2]; yel[1]=ynodes[2];
499+ edge2=SegmentIntersect(&beta1,&beta2, xel,yel,xsegment,ysegment);
500+
501+ /*edge 3: */
502+ xel[0]=xnodes[2]; yel[0]=ynodes[2]; xel[1]=xnodes[0]; yel[1]=ynodes[0];
503+ edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
504+
505+ /*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): */
506+
507+ if( (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum) ){
508+ /*This case is impossible: */
509+ _error_("error: a line cannot go through 3 different vertices!");
510+ }
511+ else if( ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum)) ){
512+
513+ /*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */
514+ if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);}
515+ if((edge2==IntersectEnum) && (edge3==IntersectEnum)) {coord1=min(beta1,gamma1); coord2=max(beta1,gamma1);}
516+ if((edge3==IntersectEnum) && (edge1==IntersectEnum)) {coord1=min(gamma1,alpha1); coord2=max(gamma1,alpha1);}
517+
518+ /*check this segment did not intersect at a vertex of the tria: */
519+ if(coord1!=coord2){
520+
521+ xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
522+ xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
523+ yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
524+ yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
525+
526+ segments_dataset->AddObject(new Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
527+ }
528+ else{
529+ /*the segment intersected at the vertex, do not bother with this "0" length segment!:*/
530+ }
531+ }
532+ else if( (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum) ){
533+
534+ /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element,
535+ * this will decide the coordinate: */
536+ if (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])){
537+ coord1=0;
538+ if(edge1==IntersectEnum){coord2=alpha1;}
539+ if(edge2==IntersectEnum){coord2=beta1;}
540+ if(edge3==IntersectEnum){coord2=gamma1;}
541+ }
542+ else{
543+ if(edge1==IntersectEnum){coord1=alpha1;}
544+ if(edge2==IntersectEnum){coord1=beta1;}
545+ if(edge3==IntersectEnum){coord1=gamma1;}
546+ coord2=1.0;
547+ }
548+
549+ xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
550+ xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
551+ yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
552+ yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
553+
554+ segments_dataset->AddObject(new Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
555+ }
556+ else{
557+ /*No interesections, but the segment might be entirely inside this triangle!: */
558+ if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){
559+ segments_dataset->AddObject(new Segment<double>(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
560+ }
561+ }
562+}/*}}}*/
563+bool NodeInElement(double* xnodes, double* ynodes, double x, double y){/*{{{*/
564+
565+ double x1,y1;
566+ double x2,y2;
567+ double x3,y3;
568+ double lambda1,lambda2,lambda3;
569+ double det;
570+
571+ x1=xnodes[0];
572+ x2=xnodes[1];
573+ x3=xnodes[2];
574+ y1=ynodes[0];
575+ y2=ynodes[1];
576+ y3=ynodes[2];
577+
578+ /*compute determinant: */
579+ det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1;
580+
581+ /*area coordinates: */
582+ lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det;
583+ lambda2=((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/det;
584+ lambda3=1-lambda1-lambda2;
585+
586+ if( ((lambda1<=1) && (lambda1>=0)) && ((lambda2<=1) && (lambda2>=0)) && ((lambda3<=1) && (lambda3>=0)) )return true;
587+ else return false;
588+
589+}/*}}}*/
590+int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){/*{{{*/
591+
592+ /*See ISSM_DIR/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */
593+
594+ /*output: */
595+ double alpha=-1;
596+ double beta=-1;
597+
598+ double xA,xB,xC,xD,yA,yB,yC,yD;
599+ double O2A[2],O2B[2],O1C[2],O1D[2];
600+ double n1[2],n2[2];
601+ double test1, test2, test3, test4;
602+ double det;
603+ double O2O1[2];
604+ double pO1A,pO1B,pO1C,pO1D;
605+
606+ xA=x1[0]; yA=y1[0];
607+ xB=x1[1]; yB=y1[1];
608+ xC=x2[0]; yC=y2[0];
609+ xD=x2[1]; yD=y2[1];
610+
611+ O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2);
612+ O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2);
613+ O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2);
614+ O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2);
615+
616+ n1[0]=yA-yB; n1[1]=xB-xA; //normal vector to segA
617+ n2[0]=yC-yD; n2[1]=xD-xC; //normal vector to segB
618+
619+ test1=n2[0]*O2A[0]+n2[1]*O2A[1];
620+ test2=n2[0]*O2B[0]+n2[1]*O2B[1];
621+
622+ if (test1*test2>0){
623+ return SeparateEnum;
624+ }
625+
626+ test3=n1[0]*O1C[0]+n1[1]*O1C[1];
627+ test4=n1[0]*O1D[0]+n1[1]*O1D[1];
628+
629+ if (test3*test4>0){
630+ return SeparateEnum;
631+ }
632+
633+ /*If colinear: */
634+ det=n1[0]*n2[1]-n2[0]*n1[1];
635+
636+ if(test1*test2==0 && test3*test4==0 && det==0){
637+
638+ //projection on the axis O1O2
639+ O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2);
640+ O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2);
641+
642+ pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]);
643+ pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]);
644+ pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1];
645+ pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1];
646+
647+ //test if one point is included in the other segment (->intersects=true)
648+ if ((pO1C-pO1A)*(pO1D-pO1A)<0){
649+ alpha=0; beta=0;
650+ *palpha=alpha;*pbeta=beta;
651+ return ColinearEnum;
652+ }
653+ if ((pO1C-pO1B)*(pO1D-pO1B)<0){
654+ alpha=0; beta=0;
655+ *palpha=alpha;*pbeta=beta;
656+ return ColinearEnum;
657+ }
658+ if ((pO1A-pO1C)*(pO1B-pO1C)<0){
659+ alpha=0; beta=0;
660+ *palpha=alpha;*pbeta=beta;
661+ return ColinearEnum;
662+ }
663+ if ((pO1A-pO1D)*(pO1B-pO1D)<0){
664+ alpha=0; beta=0;
665+ *palpha=alpha;*pbeta=beta;
666+ return ColinearEnum;
667+ }
668+
669+ //test if the 2 segments have the same middle (->intersects=true)
670+ if (O2O1==0){
671+ alpha=0; beta=0;
672+ *palpha=alpha;*pbeta=beta;
673+ return ColinearEnum;
674+ }
675+
676+ //if we are here, both segments are colinear, but do not interset:
677+ alpha=-1; beta=-1;
678+ *palpha=alpha;*pbeta=beta;
679+ return SeparateEnum;
680+ }
681+
682+ /*if we are here, both segments intersect. Determine where in the segment coordinate
683+ * system: */
684+ beta=-1;
685+ 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
686+
687+ *palpha=alpha;*pbeta=beta;
688+ return IntersectEnum;
689+}/*}}}*/
690Index: ../trunk-jpl/src/c/Makefile.am
691===================================================================
692--- ../trunk-jpl/src/c/Makefile.am (revision 14217)
693+++ ../trunk-jpl/src/c/Makefile.am (revision 14218)
694@@ -780,11 +780,6 @@
695 ./modules/AverageFilterx/AverageFilterx.h\
696 ./modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp\
697 ./modules/MeshProfileIntersectionx/MeshProfileIntersectionx.h\
698- ./modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp\
699- ./modules/MeshProfileIntersectionx/ElementSegmentsIntersection.cpp\
700- ./modules/MeshProfileIntersectionx/ElementSegment.cpp\
701- ./modules/MeshProfileIntersectionx/SegmentIntersect.cpp\
702- ./modules/MeshProfileIntersectionx/NodeInElement.cpp\
703 ./modules/ContourToMeshx/ContourToMeshx.cpp\
704 ./modules/ContourToMeshx/ContourToMeshxt.cpp\
705 ./modules/ContourToMeshx/ContourToMeshx.h\
Note: See TracBrowser for help on using the repository browser.