Ice Sheet System Model  4.18
Code documentation
MeshProfileIntersectionx.cpp
Go to the documentation of this file.
1 
4 #include "./MeshProfileIntersectionx.h"
5 
6 void MeshProfileIntersectionx(double** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods, Contour<IssmPDouble>** contours,int numcontours){/*{{{*/
7 
8  int i,j,k;
9 
10  /*Contour:*/
11  Contour<IssmPDouble>* contouri=NULL;
12  int numnodes;
13  double* xc=NULL;
14  double* yc=NULL;
15 
16  /*output: */
17  double* segments=NULL;
18  int numsegs;
19 
20  /*intermediary: */
21  double** allsegments=NULL;
22  double* segmentsi=NULL;
23  int* allnumsegs=NULL;
24  int numsegsi;
25  int count;
26 
27  /*Allocate: */
28  allsegments=xNew<double*>(numcontours);
29  allnumsegs=xNew<int>(numcontours);
30 
31  /*Loop through all contours: */
32  for (i=0;i<numcontours;i++){
33 
34  /*retrieve contour info: */
35  contouri=*(contours+i);
36  numnodes=contouri->nods;
37  xc=contouri->x;
38  yc=contouri->y;
39 
40  /*determine segmentsi and numsegsi for this contour and the mesh intersection: */
41  MeshSegmentsIntersection(&segmentsi,&numsegsi,index,x,y,nel,nods,xc,yc,numnodes);
42 
43  /*save segmentsi: */
44  allsegments[i]=segmentsi;
45  allnumsegs[i]=numsegsi;
46  }
47 
48  /*total number of segments? */
49  numsegs=0;
50  for(i=0;i<numcontours;i++)numsegs+=allnumsegs[i];
51 
52  /*Out of all segments, create one common array of segments: */
53  segments=xNew<double>(5*numsegs);
54  count=0;
55  for(i=0;i<numcontours;i++){
56 
57  segmentsi=allsegments[i];
58  numsegsi=allnumsegs[i];
59 
60  for(j=0;j<numsegsi;j++){
61  for(k=0;k<5;k++){
62  *(segments+count*5+k)=*(segmentsi+j*5+k);
63  }
64  count++;
65  }
66  }
67 
68  /*Assign output pointers:*/
69  *psegments=segments;
70  *pnumsegs=numsegs;
71 }/*}}}*/
72 void 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  int numsegs;
79 
80  /*intermediary: */
81  DataSet* segments_dataset=NULL;
82  double xnodes[3];
83  double ynodes[3];
84 
85  /*We don't know how many segments we are going to get, so have a dynamic container: */
86  segments_dataset=new DataSet();
87 
88  /*Go through elements, and call ElementSegmentsIntersection routine: */
89  for(i=0;i<nel;i++){
90  for(j=0;j<3;j++){
91  xnodes[j]=x[*(index+3*i+j)];
92  ynodes[j]=y[*(index+3*i+j)];
93  }
94  ElementSegmentsIntersection(segments_dataset,i,xnodes,ynodes,xc,yc,numnodes);
95  }
96 
97  /*Using the segments_dataset dataset, create segments: */
98  numsegs=segments_dataset->Size();
99  segments=xNew<double>(5*numsegs);
100  for(i=0;i<numsegs;i++){
101  Segment<double>* segment=(Segment<double>*)segments_dataset->GetObjectByOffset(i);
102 
103  /*x1,y1,x2,y2 then element_id: */
104  segments[5*i+0]=segment->x1;
105  segments[5*i+1]=segment->y1;
106  segments[5*i+2]=segment->x2;
107  segments[5*i+3]=segment->y2;
108  segments[5*i+4]=(double)segment->eid;
109  }
110 
111  /*Free ressources:*/
112  delete segments_dataset;
113 
114  /*Assign output pointers:*/
115  *psegments=segments;
116  *pnumsegs=numsegs;
117 }/*}}}*/
118 
119 /*Utilities*/
120 void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes){/*{{{*/
121 
122  double xsegment[2];
123  double ysegment[2];
124 
125  /*Loop through contour: */
126  for(int i=0;i<numnodes-1;i++){
127  xsegment[0]=xc[i];
128  xsegment[1]=xc[i+1];
129  ysegment[0]=yc[i];
130  ysegment[1]=yc[i+1];
131  /*if (el==318 && i==9){
132  _printf_("contour: " << i << " " << xsegment[0] << " " << ysegment[0] << " " << xsegment[1] << " " << ysegment[1]
133  << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " <<
134  ynodes[2] << "\n");
135  }*/
136  ElementSegment(segments_dataset,el, i, xnodes,ynodes,xsegment,ysegment);
137  }
138 }/*}}}*/
139 void ElementSegment(DataSet* segments_dataset,int el, int contouri, double* xnodes,double* ynodes,double* xsegment,double* ysegment){/*{{{*/
140 
141  /*We have a tria element (xnodes,ynodes) and a segment (xsegment,ysegment). Find whether they intersect.
142  * If they do, create a Segment object with the intersection, and add to segments_dataset dataset: */
143 
144  double alpha1,alpha2;
145  double beta1,beta2;
146  double gamma1,gamma2;
147 
148  int edge1,edge2,edge3;
149 
150  double xel[2],yel[2];
151  double coord1 = 0.;
152  double coord2 = 0.;
153  double xfinal[2],yfinal[2];
154 
155  /*edge 1: */
156  xel[0]=xnodes[0]; yel[0]=ynodes[0]; xel[1]=xnodes[1]; yel[1]=ynodes[1];
157  edge1=SegmentIntersect(&alpha1,&alpha2, xel,yel,xsegment,ysegment); //alpha1: segment coordinate of intersection. alpha2: same thing for second interesection if it exists (colinear edges)
158 
159  /*edge 2: */
160  xel[0]=xnodes[1]; yel[0]=ynodes[1]; xel[1]=xnodes[2]; yel[1]=ynodes[2];
161  edge2=SegmentIntersect(&beta1,&beta2, xel,yel,xsegment,ysegment);
162 
163  /*edge 3: */
164  xel[0]=xnodes[2]; yel[0]=ynodes[2]; xel[1]=xnodes[0]; yel[1]=ynodes[0];
165  edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
166 
167  /*edge can be either IntersectEnum (one and only one intersection between the edge and the segment), ColinearEnum (edge and segment are collinear) and SeparateEnum (no intersection): */
168 
169  /*if (el==318 && contouri==9){
170  _printf_(edge1 << " " << edge2 << " " << edge3 << " " << alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " " << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
171 
172  _printf_("Bool" << (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum));
173  }*/
174 
175  if( (edge1==IntersectEnum) && (edge2==IntersectEnum) && (edge3==IntersectEnum) ){
176 
177  /*This can only be the case if the segment intersected through one vertex, meaning a pair from alpha1, beta1 or gamma1 is 0:*/
178  if (alpha1!=0 && alpha1!=1){
179  /*The vertex opposite edge 1 was intersected:*/
180  xfinal[0]=xsegment[0]+alpha1*(xsegment[1]-xsegment[0]);
181  yfinal[0]=ysegment[0]+alpha1*(ysegment[1]-ysegment[0]);
182  xfinal[1]=xnodes[2];
183  yfinal[1]=ynodes[2];
184  }
185  else if (beta1!=0 && beta1!=1){
186  /*The vertex opposite edge 2 was intersected:*/
187  xfinal[0]=xsegment[0]+beta1*(xsegment[1]-xsegment[0]);
188  yfinal[0]=ysegment[0]+beta1*(ysegment[1]-ysegment[0]);
189  xfinal[1]=xnodes[0];
190  yfinal[1]=ynodes[0];
191  }
192  else if (gamma1!=0 && gamma1!=1){
193  /*The vertex opposite edge 3 was intersected:*/
194  xfinal[0]=xsegment[0]+gamma1*(xsegment[1]-xsegment[0]);
195  yfinal[0]=ysegment[0]+gamma1*(ysegment[1]-ysegment[0]);
196  xfinal[1]=xnodes[1];
197  yfinal[1]=ynodes[1];
198  }
199  segments_dataset->AddObject(new Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
200 
201  /*This case is impossible: not quite! */
202  //_printf_(alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " " << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
203  /* _error_("error: a line cannot go through 3 different vertices!");*/
204  }
205  else if( ((edge1==IntersectEnum) && (edge2==IntersectEnum)) || ((edge2==IntersectEnum) && (edge3==IntersectEnum)) || ((edge3==IntersectEnum) && (edge1==IntersectEnum)) ){
206 
207  /*segment interscts 2 opposite edges of our triangle, at 2 segment coordinates, pick up the lowest (coord1) and highest (coord2): */
208  if((edge1==IntersectEnum) && (edge2==IntersectEnum)) {coord1=min(alpha1,beta1); coord2=max(alpha1,beta1);}
209  if((edge2==IntersectEnum) && (edge3==IntersectEnum)) {coord1=min(beta1,gamma1); coord2=max(beta1,gamma1);}
210  if((edge3==IntersectEnum) && (edge1==IntersectEnum)) {coord1=min(gamma1,alpha1); coord2=max(gamma1,alpha1);}
211 
212  /*check this segment did not intersect at a vertex of the tria: */
213  if(coord1!=coord2){
214 
215  xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
216  xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
217  yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
218  yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
219 
220  segments_dataset->AddObject(new Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
221  }
222  else{
223  /*the segment intersected at the vertex, do not bother with this "0" length segment!:*/
224  }
225  }
226  else if( (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum) ){
227 
228  /*if (el==318 && contouri==9){
229  _printf_("hello" << " NodeInElement 0 " << (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) << " NodeInElement 1 " << (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])));
230  }*/
231 
232  /*segment intersect only 1 edge. Figure out where the first point in the segment is, inside or outside the element,
233  * this will decide the coordinate: */
234  if (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])){
235  coord1=0;
236  if(edge1==IntersectEnum){coord2=alpha1;}
237  if(edge2==IntersectEnum){coord2=beta1;}
238  if(edge3==IntersectEnum){coord2=gamma1;}
239  }
240  else if (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])){
241  if(edge1==IntersectEnum){coord1=alpha1;}
242  if(edge2==IntersectEnum){coord1=beta1;}
243  if(edge3==IntersectEnum){coord1=gamma1;}
244  coord2=1.0;
245  }
246  else{
247  double tolerance=1e-10;
248  /*Ok, we have an issue here. Probably one of the segments' end is on a vertex, within a certain tolerance!*/
249  if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[0],ysegment[0],tolerance) ||
250  IsIdenticalNode(xnodes[1],ynodes[1],xsegment[0],ysegment[0],tolerance) ||
251  IsIdenticalNode(xnodes[2],ynodes[2],xsegment[0],ysegment[0],tolerance)){
252 
253  /*ok, segments[0] is common to one of our vertices: */
254  coord1=0;
255  if(edge1==IntersectEnum){coord2=alpha1;}
256  if(edge2==IntersectEnum){coord2=beta1;}
257  if(edge3==IntersectEnum){coord2=gamma1;}
258  }
259  else if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[1],ysegment[1],tolerance) ||
260  IsIdenticalNode(xnodes[1],ynodes[1],xsegment[1],ysegment[1],tolerance) ||
261  IsIdenticalNode(xnodes[2],ynodes[2],xsegment[1],ysegment[1],tolerance)){
262 
263  /*ok, segments[1] is common to one of our vertices: */
264  //if (el==318 && contouri==9){ _printf_("ok2" << "\n"); }
265  if(edge1==IntersectEnum){coord1=alpha1;}
266  if(edge2==IntersectEnum){coord1=beta1;}
267  if(edge3==IntersectEnum){coord1=gamma1;}
268  coord2=1.0;
269  }
270  }
271 
272  xfinal[0]=xsegment[0]+coord1*(xsegment[1]-xsegment[0]);
273  xfinal[1]=xsegment[0]+coord2*(xsegment[1]-xsegment[0]);
274  yfinal[0]=ysegment[0]+coord1*(ysegment[1]-ysegment[0]);
275  yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
276 
277  segments_dataset->AddObject(new Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
278  }
279  else{
280  /*No interesections, but the segment might be entirely inside this triangle!: */
281  if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){
282  segments_dataset->AddObject(new Segment<double>(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
283  }
284  }
285 }/*}}}*/
286 bool NodeInElement(double* xnodes, double* ynodes, double x, double y){/*{{{*/
287 
288  double x1,y1;
289  double x2,y2;
290  double x3,y3;
291  double lambda1,lambda2,lambda3;
292  double det;
293 
294  x1=xnodes[0];
295  x2=xnodes[1];
296  x3=xnodes[2];
297  y1=ynodes[0];
298  y2=ynodes[1];
299  y3=ynodes[2];
300 
301  /*compute determinant: */
302  det=x1*y2-x1*y3-x3*y2-x2*y1+x2*y3+x3*y1;
303 
304  /*area coordinates: */
305  lambda1=((y2-y3)*(x-x3)+(x3-x2)*(y-y3))/det;
306  lambda2=((y3-y1)*(x-x3)+(x1-x3)*(y-y3))/det;
307  lambda3=1-lambda1-lambda2;
308 
309  if( ((lambda1<=1) && (lambda1>=0)) && ((lambda2<=1) && (lambda2>=0)) && ((lambda3<=1) && (lambda3>=0)) )return true;
310  else return false;
311 
312 }/*}}}*/
313 int SegmentIntersect(double* palpha, double* pbeta, double* x1, double* y1, double* x2, double* y2){/*{{{*/
314 
315  /*See ISSM_DIR/src/m/utils/Geometry/SegIntersect.m for matlab routine from which we take this routine: */
316 
317  /*output: */
318  double alpha=-1;
319  double beta=-1;
320 
321  double xA,xB,xC,xD,yA,yB,yC,yD;
322  double O2A[2],O2B[2],O1C[2],O1D[2];
323  double n1[2],n2[2];
324  double test1, test2, test3, test4;
325  double det;
326  double O2O1[2];
327  double pO1A,pO1B,pO1C,pO1D;
328 
329  xA=x1[0]; yA=y1[0];
330  xB=x1[1]; yB=y1[1];
331  xC=x2[0]; yC=y2[0];
332  xD=x2[1]; yD=y2[1];
333 
334  O2A[0]=xA -(xD/2+xC/2); O2A[1]=yA -(yD/2+yC/2);
335  O2B[0]=xB -(xD/2+xC/2); O2B[1]=yB -(yD/2+yC/2);
336  O1C[0]=xC -(xA/2+xB/2); O1C[1]=yC -(yA/2+yB/2);
337  O1D[0]=xD -(xA/2+xB/2); O1D[1]=yD -(yA/2+yB/2);
338 
339  n1[0]=yA-yB; n1[1]=xB-xA; //normal vector to segA
340  n2[0]=yC-yD; n2[1]=xD-xC; //normal vector to segB
341 
342  test1=n2[0]*O2A[0]+n2[1]*O2A[1];
343  test2=n2[0]*O2B[0]+n2[1]*O2B[1];
344 
345  if (test1*test2>0){
346  return SeparateEnum;
347  }
348 
349  test3=n1[0]*O1C[0]+n1[1]*O1C[1];
350  test4=n1[0]*O1D[0]+n1[1]*O1D[1];
351 
352  if (test3*test4>0){
353  return SeparateEnum;
354  }
355 
356  /*If colinear: */
357  det=n1[0]*n2[1]-n2[0]*n1[1];
358 
359  if(test1*test2==0 && test3*test4==0 && det==0){
360 
361  //projection on the axis O1O2
362  O2O1[0]=(xA/2+xB/2)-(xD/2+xC/2);
363  O2O1[1]=(yA/2+yB/2)-(yD/2+yC/2);
364 
365  pO1A=O2O1[0]*(O2A[0]-O2O1[0])+O2O1[1]*(O2A[1]-O2O1[1]);
366  pO1B=O2O1[0]*(O2B[0]-O2O1[0])+O2O1[1]*(O2B[1]-O2O1[1]);
367  pO1C=O2O1[0]*O1C[0]+O2O1[1]*O1C[1];
368  pO1D=O2O1[0]*O1D[0]+O2O1[1]*O1D[1];
369 
370  //test if one point is included in the other segment (->intersects=true)
371  if ((pO1C-pO1A)*(pO1D-pO1A)<0){
372  alpha=0; beta=0;
373  *palpha=alpha;*pbeta=beta;
374  return ColinearEnum;
375  }
376  if ((pO1C-pO1B)*(pO1D-pO1B)<0){
377  alpha=0; beta=0;
378  *palpha=alpha;*pbeta=beta;
379  return ColinearEnum;
380  }
381  if ((pO1A-pO1C)*(pO1B-pO1C)<0){
382  alpha=0; beta=0;
383  *palpha=alpha;*pbeta=beta;
384  return ColinearEnum;
385  }
386  if ((pO1A-pO1D)*(pO1B-pO1D)<0){
387  alpha=0; beta=0;
388  *palpha=alpha;*pbeta=beta;
389  return ColinearEnum;
390  }
391 
392  //test if the 2 segments have the same middle (->intersects=true)
393  if(O2O1[0]==0 && O2O1[1]){
394  alpha=0; beta=0;
395  *palpha=alpha;*pbeta=beta;
396  return ColinearEnum;
397  }
398 
399  //if we are here, both segments are colinear, but do not interset:
400  alpha=-1; beta=-1;
401  *palpha=alpha;*pbeta=beta;
402  return SeparateEnum;
403  }
404 
405  /*if we are here, both segments intersect. Determine where in the segment coordinate
406  * system: */
407  beta=-1;
408  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
409 
410  *palpha=alpha;*pbeta=beta;
411  return IntersectEnum;
412 }/*}}}*/
413 bool IsIdenticalNode(double x1, double y1, double x2, double y2, double tolerance){ /*{{{*/
414 
415  if (sqrt(pow(x1-x2,2.0) + pow(y1-y2,2))<tolerance)return true;
416  else return false;
417 
418 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
SegmentIntersect
int SegmentIntersect(double *palpha, double *pbeta, double *x1, double *y1, double *x2, double *y2)
Definition: MeshProfileIntersectionx.cpp:313
IntersectEnum
@ IntersectEnum
Definition: EnumDefinitions.h:1132
SeparateEnum
@ SeparateEnum
Definition: EnumDefinitions.h:1272
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
Segment::x1
doubletype x1
Definition: Segment.h:19
Contour
Definition: Contour.h:15
Segment
Definition: Segment.h:15
NodeInElement
bool NodeInElement(double *xnodes, double *ynodes, double x, double y)
Definition: MeshProfileIntersectionx.cpp:286
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
bamg::det
long long det(const I2 &a, const I2 &b, const I2 &c)
Definition: det.h:8
MeshSegmentsIntersection
void MeshSegmentsIntersection(double **psegments, int *pnumsegs, int *index, double *x, double *y, int nel, int nods, double *xc, double *yc, int numnodes)
Definition: MeshProfileIntersectionx.cpp:72
Segment::eid
int eid
Definition: Segment.h:18
Segment::y1
doubletype y1
Definition: Segment.h:20
Contour::nods
int nods
Definition: Contour.h:20
Segment::y2
doubletype y2
Definition: Segment.h:22
MeshProfileIntersectionx
void MeshProfileIntersectionx(double **psegments, int *pnumsegs, int *index, double *x, double *y, int nel, int nods, Contour< IssmPDouble > **contours, int numcontours)
Definition: MeshProfileIntersectionx.cpp:6
Contour::y
doubletype * y
Definition: Contour.h:22
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
ColinearEnum
@ ColinearEnum
Definition: EnumDefinitions.h:1012
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
Contour::x
doubletype * x
Definition: Contour.h:21
ElementSegmentsIntersection
void ElementSegmentsIntersection(DataSet *segments_dataset, int el, double *xnodes, double *ynodes, double *xc, double *yc, int numnodes)
Definition: MeshProfileIntersectionx.cpp:120
Segment::x2
doubletype x2
Definition: Segment.h:21
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
MeshProfileIntersectionx.h
DataSet
Declaration of DataSet class.
Definition: DataSet.h:14
ElementSegment
void ElementSegment(DataSet *segments_dataset, int el, int contouri, double *xnodes, double *ynodes, double *xsegment, double *ysegment)
Definition: MeshProfileIntersectionx.cpp:139
IsIdenticalNode
bool IsIdenticalNode(double x1, double y1, double x2, double y2, double tolerance)
Definition: MeshProfileIntersectionx.cpp:413
Generated on Thu Jul 2 2020 08:09:19 for Ice Sheet System Model by   doxygen 1.8.19