Ice Sheet System Model  4.18
Code documentation
ListofIntersectionTriangles.cpp
Go to the documentation of this file.
1 #include <cstdio>
2 #include <cstring>
3 #include <cmath>
4 #include <ctime>
5 
6 #include "./bamgobjects.h"
7 #include "../shared/shared.h"
8 #include "./det.h"
9 
10 namespace bamg {
11 
12  /*Constructors Destructors*/
13  ListofIntersectionTriangles::ListofIntersectionTriangles(int n,int m)/*{{{*/
14  : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
15  NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){
16  }
17  /*}}}*/
18  ListofIntersectionTriangles::~ListofIntersectionTriangles(){/*{{{*/
19  if (lIntTria) delete [] lIntTria,lIntTria=0;
20  if (lSegsI) delete [] lSegsI,lSegsI=0;
21  }
22  /*}}}*/
23 
24  /*Methods*/
25  void ListofIntersectionTriangles::Init(void){/*{{{*/
26  state=0;
27  len=0;
28  Size=0;
29  }
30  /*}}}*/
31  double ListofIntersectionTriangles::Length(){/*{{{*/
32  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Length)*/
33 
34  // computation of the length
35 
36  // check Size
37  if (Size<=0){
38  _error_("Size<=0");
39  }
40 
41  Metric Mx,My;
42  int ii,jj;
43  R2 x,y,xy;
44 
45  SegInterpolation* SegI=lSegsI;
46  lSegsI[NbSeg].last=Size;
47  int EndSeg=Size;
48 
49  y = lIntTria[0].x;
50  double sxy, s = 0;
51  lIntTria[0].s =0;
52  SegI->lBegin=s;
53 
54  for (jj=0,ii=1;ii<Size;jj=ii++) {
55  // seg jj,ii
56  x = y;
57  y = lIntTria[ii].x;
58  xy = y-x;
59  Mx = lIntTria[ii].m;
60  My = lIntTria[jj].m;
61  sxy= LengthInterpole(Mx,My,xy);
62  s += sxy;
63  lIntTria[ii].s = s;
64  if (ii == EndSeg){
65  SegI->lEnd=s;
66  SegI++;
67  EndSeg=SegI->last;
68  SegI->lBegin=s;
69  }
70  }
71  len = s;
72  SegI->lEnd=s;
73 
74  return s;
75  }
76  /*}}}*/
77  int ListofIntersectionTriangles::NewItem(Triangle * tt,double d0,double d1,double d2) { /*{{{*/
78  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
79 
80  int n;
81  R2 x(0,0);
82  if ( d0) x = (*tt)[0].r * d0;
83  if ( d1) x = x + (*tt)[1].r * d1;
84  if ( d2) x = x + (*tt)[2].r * d2;
85  // newer add same point
86  if(!Size || Norme2_2(lIntTria[Size-1].x-x)) {
87  if (Size==MaxSize) ReShape();
88  lIntTria[Size].t=tt;
89  lIntTria[Size].bary[0]=d0;
90  lIntTria[Size].bary[1]=d1;
91  lIntTria[Size].bary[2]=d2;
92  lIntTria[Size].x = x;
93  Metric m0,m1,m2;
94  BamgVertex * v;
95  if ((v=(*tt)(0))) m0 = v->m;
96  if ((v=(*tt)(1))) m1 = v->m;
97  if ((v=(*tt)(2))) m2 = v->m;
98  lIntTria[Size].m = Metric(lIntTria[Size].bary,m0,m1,m2);
99  n=Size++;}
100  else n=Size-1;
101  return n;
102  }
103  /*}}}*/
104  int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {/*{{{*/
105  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewItem)*/
106 
107  int n;
108  if(!Size || Norme2_2(lIntTria[Size-1].x-A)) {
109  if (Size==MaxSize) ReShape();
110  lIntTria[Size].t=0;
111  lIntTria[Size].x=A;
112  lIntTria[Size].m=mm;
113  n=Size++;
114  }
115  else n=Size-1;
116  return n;
117  }
118  /*}}}*/
119  long ListofIntersectionTriangles::NewPoints(BamgVertex* vertices,long &nbv,long maxnbv){/*{{{*/
120  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/NewPoints)*/
121 
122  //If length<1.5, do nothing
123  double s=Length();
124  if (s<1.5) return 0;
125 
126  const long nbvold=nbv;
127  int ii = 1 ;
128  R2 y,x;
129  Metric My,Mx ;
130  double sx =0,sy;
131  int nbi=Max(2,(int) (s+0.5));
132  double sint=s/nbi;
133  double si =sint;
134 
135  int EndSeg=Size;
136  SegInterpolation* SegI=NULL;
137  if (NbSeg) SegI=lSegsI,EndSeg=SegI->last;
138 
139  for (int k=1;k<nbi;k++){
140  while ((ii < Size) && ( lIntTria[ii].s <= si )){
141  if (ii++ == EndSeg){
142  SegI++;
143  EndSeg=SegI->last;
144  }
145  }
146 
147  int ii1=ii-1;
148  x =lIntTria[ii1].x;
149  sx =lIntTria[ii1].s;
150  Metric Mx=lIntTria[ii1].m;
151  y =lIntTria[ii].x;
152  sy =lIntTria[ii].s;
153  Metric My=lIntTria[ii].m;
154  double lxy = sy-sx;
155  double cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);
156 
157  R2 C;
158  double cx = 1-cy;
159  C = SegI ? SegI->F(si): x * cx + y *cy;
160  //C.Echo();
161  //x.Echo();
162  //y.Echo();
163  //_printf_("cx = " << cx << ", cy=" << cy << "\n");
164 
165  si += sint;
166  if ( nbv<maxnbv) {
167  vertices[nbv].r = C;
168  vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);
169  }
170  else return nbv-nbvold;
171  }
172  return nbv-nbvold;
173  }
174  /*}}}*/
175  void ListofIntersectionTriangles::ReShape(){ /*{{{*/
176 
177  int newsize = MaxSize*2;
178  IntersectionTriangles* nw = new IntersectionTriangles[newsize];
179  _assert_(nw);
180 
181  // recopy
182  for (int i=0;i<MaxSize;i++) nw[i] = lIntTria[i];
183  long int verbosity=0;
184  if(verbosity>3) _printf_(" ListofIntersectionTriangles ReShape Maxsize " << MaxSize << " -> " << MaxNbSeg << "\n");
185  MaxSize = newsize;
186  delete [] lIntTria;// remove old
187  lIntTria = nw; // copy pointer
188  }
189  /*}}}*/
190  void ListofIntersectionTriangles::SplitEdge(Mesh & Bh, const R2 &A,const R2 &B,int nbegin) {/*{{{*/
191  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ListofIntersectionTriangles)*/
192 
193  Triangle *tbegin, *t;
194 
195  long long deta[3], deti,detj;
196  double ba[3];
197  int ifirst=-1,ilast;
198  int i0,i1,i2;
199  int ocut,i,j,k=-1;
200  // int OnAVertices =0;
201  long long dt[3];
202  I2 a = Bh.R2ToI2(A) ,b= Bh.R2ToI2(B);// compute the Icoor a,b
203  I2 vi,vj;
204  int iedge =-1;// not a edge
205 
206  if(nbegin) {// optimisation
207  // we suppose knowing the starting triangle
208  t=tbegin=lIntTria[ilast=(Size-1)].t;
209  if (tbegin->det>=0)
210  ifirst = ilast;}
211  else {// not optimisation
212  Init();
213  t=tbegin = Bh.TriangleFindFromCoord(a,deta);
214  if( t->det>=0)
215  ilast=NewItem(t,double(deta[0])/t->det,double(deta[1])/t->det,double(deta[2])/t->det);
216  else
217  {// find the nearest boundary edge of the vertex A
218  // find a edge or such normal projection a the edge IJ is on the edge
219  // <=> IJ.IA >=0 && IJ.AJ >=0
220  ilast=ifirst;
221  double ba,bb;
222  AdjacentTriangle edge=CloseBoundaryEdge(a,t,ba,bb);
223  BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
224  NewItem(A,Metric(ba,v0,bb,v1));
225  t=edge;
226  // test if the point b is in the same side
227  if (det(v0.i,v1.i,b)>=0) {
228  AdjacentTriangle edge=CloseBoundaryEdge(a,t,ba,bb);
229  BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
230  NewItem(A,Metric(ba,v0,bb,v1));
231  return;
232  }
233  } // find the nearest boundary edge of the vertex A
234  } // end not optimisation
235  if (t->det<0) { // outside departure
236  while (t->det <0) { // intersection boundary edge and a,b,
237  k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;
238  if (k<0){
239  _error_("k<0");
240  }
241  ocut = OppositeEdge[k];
242  i=VerticesOfTriangularEdge[ocut][0];
243  j=VerticesOfTriangularEdge[ocut][1];
244  vi=(*t)[i];
245  vj=(*t)[j];
246  deti = bamg::det(a,b,vi);
247  detj = bamg::det(a,b,vj);
248  if (deti>0) // go to i direction on gamma
249  ocut = PreviousEdge[ocut];
250  else if (detj<=0) // go to j direction on gamma
251  ocut = NextEdge[ocut];
252  AdjacentTriangle tadj =t->Adj(ocut);
253  t = tadj;
254  iedge= tadj;
255  if (t == tbegin) { //
256  double ba,bb;
257  AdjacentTriangle edge=CloseBoundaryEdge(a,t,ba,bb);
258  BamgVertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);
259  NewItem(A,Metric(ba,v0,bb,v1));
260  return;
261  }
262  } // end while (t->det <0)
263  // theoriticaly we have: deti =<0 and detj>0
264 
265  // computation of barycentric coor
266  // test if the point b is on size on t
267  // we revert vi,vj because vi,vj is def in Adj triangle
268  if ( det(vi,vj,b)>=0) {
269  t=tbegin;
270  double ba,bb;
271  AdjacentTriangle edge=CloseBoundaryEdge(b,t,ba,bb);
272  NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
273  return;
274  }
275  else
276  {
277  k = OppositeVertex[iedge];
278  i=VerticesOfTriangularEdge[iedge][0];
279  j=VerticesOfTriangularEdge[iedge][1];
280  double dij = detj-deti;
281  if (i+j+k != 0 + 1 +2){
282  _error_("i+j+k != 0 + 1 +2");
283  }
284  ba[j] = detj/dij;
285  ba[i] = -deti/dij;
286  ba[k] = 0;
287  ilast=NewItem(t,ba[0],ba[1],ba[2]); }
288  } // outside departure
289 
290  // recherche the intersection of [a,b] with Bh Mesh.
291  // we know a triangle ta contening the vertex a
292  // we have 2 case for intersection [a,b] with a edge [A,B] of Bh
293  // 1) the intersection point is in ]A,B[
294  // 2) is A or B
295  // first version ---
296  for (;;) {
297  // t->Draw();
298  if (iedge < 0) {
299  i0 =0;i1=1;i2=2;
300  dt[0] =bamg::det(a,b,(*t)[0]);
301  dt[1] =bamg::det(a,b,(*t)[1]);
302  dt[2] =bamg::det(a,b,(*t)[2]);}
303  else {
304  i2 = iedge;
305  i0 = NextEdge[i2];
306  i1 = NextEdge[i0];
307  dt[VerticesOfTriangularEdge[iedge][0]] = detj;// we revert i,j because
308  dt[VerticesOfTriangularEdge[iedge][1]] = deti;// we take the Triangle by the other side
309  dt[iedge] = det(a,b,(*t)[OppositeVertex[iedge]]);}
310 
311  // so we have just to see the transition from - to + of the det0..2 on edge of t
312  // because we are going from a to b
313  if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&
314  ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0))
315  ocut =i0;
316  else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&
317  (dt[j=VerticesOfTriangularEdge[i1][1]] > 0))
318  ocut =i1;
319  else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
320  (dt[j=VerticesOfTriangularEdge[i2][1]] > 0))
321  ocut =i2;
322  else if ((dt[i=VerticesOfTriangularEdge[i0][0]] == 0) &&
323  ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0))
324  ocut =i0;
325  else if ((dt[i=VerticesOfTriangularEdge[i1][0]] == 0) &&
326  (dt[j=VerticesOfTriangularEdge[i1][1]] > 0))
327  ocut =i1;
328  else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) &&
329  (dt[j=VerticesOfTriangularEdge[i2][1]] > 0))
330  ocut =i2;
331  else if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&
332  ( dt[j=VerticesOfTriangularEdge[i0][1]] == 0))
333  ocut =i0;
334  else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&
335  (dt[j=VerticesOfTriangularEdge[i1][1]] == 0))
336  ocut =i1;
337  else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&
338  (dt[j=VerticesOfTriangularEdge[i2][1]] == 0))
339  ocut =i2;
340  else { // On a edge (2 zero)
341  k =0;
342  if (dt[0]) ocut=0,k++;
343  if (dt[1]) ocut=1,k++;
344  if (dt[2]) ocut=2,k++;
345  if(k == 1) {
346  if (dt[ocut] >0) // triangle upper AB
347  ocut = NextEdge[ocut];
348  i= VerticesOfTriangularEdge[ocut][0];
349  j= VerticesOfTriangularEdge[ocut][1];
350  }
351  else {
352  _error_("Bug Split Edge");
353  }
354  }
355 
356  k = OppositeVertex[ocut];
357 
358  long long detbij = bamg::det((*t)[i],(*t)[j],b);
359 
360  if (detbij >= 0) { //we find the triangle contening b
361  dt[0]=bamg::det((*t)[1],(*t)[2],b);
362  dt[1]=bamg::det((*t)[2],(*t)[0],b);
363  dt[2]=bamg::det((*t)[0],(*t)[1],b);
364  double dd = t->det;
365  NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);
366  return ;}
367  else { // next triangle by adjacent by edge ocut
368  deti = dt[i];
369  detj = dt[j];
370  double dij = detj-deti;
371  ba[i] = detj/dij;
372  ba[j] = -deti/dij;
373  ba[3-i-j ] = 0;
374  ilast=NewItem(t, ba[0],ba[1],ba[2]);
375 
376  AdjacentTriangle ta =t->Adj(ocut);
377  t = ta;
378  iedge= ta;
379  if (t->det <= 0) {
380  double ba,bb;
381  AdjacentTriangle edge=CloseBoundaryEdge(b,t,ba,bb);
382  NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));
383  return;
384  }
385  }// we go outside of omega
386  } // for(;;)
387  }
388  /*}}}*/
389 
390 }
bamg::BamgVertex
Definition: BamgVertex.h:15
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
bamg::BamgVertex::i
I2 i
Definition: BamgVertex.h:20
bamg::BamgVertex::r
R2 r
Definition: BamgVertex.h:21
bamg::PreviousEdge
static const short PreviousEdge[3]
Definition: macros.h:18
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
bamg
Definition: AdjacentTriangle.cpp:9
bamg::abscisseInterpole
double abscisseInterpole(const Metric &Ma, const Metric &Mb, R2 AB, double s, int optim)
Definition: Metric.cpp:327
bamg::Triangle
Definition: Triangle.h:13
bamg::ListofIntersectionTriangles::IntersectionTriangles::bary
double bary[3]
Definition: ListofIntersectionTriangles.h:16
bamg::Max
T Max(const T &a, const T &b)
Definition: extrema.h:7
bamg::ListofIntersectionTriangles::IntersectionTriangles::x
R2 x
Definition: ListofIntersectionTriangles.h:17
bamg::AdjacentTriangle
Definition: AdjacentTriangle.h:12
bamg::ListofIntersectionTriangles::IntersectionTriangles::s
double s
Definition: ListofIntersectionTriangles.h:19
bamg::OppositeEdge
static const short OppositeEdge[3]
Definition: macros.h:16
bamg::Triangle::Adj
AdjacentTriangle Adj(int i) const
Definition: Triangle.cpp:48
bamg::ListofIntersectionTriangles::lIntTria
IntersectionTriangles * lIntTria
Definition: ListofIntersectionTriangles.h:48
bamg::ListofIntersectionTriangles::NewItem
int NewItem(Triangle *tt, double d0, double d1, double d2)
Definition: ListofIntersectionTriangles.cpp:77
bamg::ListofIntersectionTriangles::SegInterpolation::lEnd
double lEnd
Definition: ListofIntersectionTriangles.h:29
det.h
bamg::det
long long det(const I2 &a, const I2 &b, const I2 &c)
Definition: det.h:8
bamg::ListofIntersectionTriangles::IntersectionTriangles
Definition: ListofIntersectionTriangles.h:12
bamg::ListofIntersectionTriangles::Init
void Init()
Definition: ListofIntersectionTriangles.cpp:25
bamg::ListofIntersectionTriangles::IntersectionTriangles::t
Triangle * t
Definition: ListofIntersectionTriangles.h:15
bamg::ListofIntersectionTriangles::NewPoints
long NewPoints(BamgVertex *, long &nbv, long maxnbv)
Definition: ListofIntersectionTriangles.cpp:119
bamg::Metric
Definition: Metric.h:17
bamg::ListofIntersectionTriangles::Size
int Size
Definition: ListofIntersectionTriangles.h:45
bamg::LengthInterpole
double LengthInterpole(const Metric &Ma, const Metric &Mb, R2 AB)
Definition: Metric.cpp:158
bamg::ListofIntersectionTriangles::~ListofIntersectionTriangles
~ListofIntersectionTriangles()
Definition: ListofIntersectionTriangles.cpp:18
bamg::ListofIntersectionTriangles::ListofIntersectionTriangles
ListofIntersectionTriangles(int n=256, int m=16)
Definition: ListofIntersectionTriangles.cpp:13
bamg::ListofIntersectionTriangles::MaxNbSeg
int MaxNbSeg
Definition: ListofIntersectionTriangles.h:50
bamg::ListofIntersectionTriangles::state
int state
Definition: ListofIntersectionTriangles.h:47
bamg::ListofIntersectionTriangles::SegInterpolation
Definition: ListofIntersectionTriangles.h:24
bamg::Mesh
Definition: Mesh.h:21
bamg::ListofIntersectionTriangles::IntersectionTriangles::m
Metric m
Definition: ListofIntersectionTriangles.h:18
bamg::ListofIntersectionTriangles::lSegsI
SegInterpolation * lSegsI
Definition: ListofIntersectionTriangles.h:51
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
bamg::Triangle::det
long long det
Definition: Triangle.h:23
bamg::CloseBoundaryEdge
AdjacentTriangle CloseBoundaryEdge(I2 A, Triangle *t, double &a, double &b)
Definition: Mesh.cpp:4927
bamg::VerticesOfTriangularEdge
static const short VerticesOfTriangularEdge[3][2]
Definition: macros.h:13
bamg::Mesh::R2ToI2
I2 R2ToI2(const R2 &P) const
Definition: Mesh.cpp:3937
bamg::ListofIntersectionTriangles::len
double len
Definition: ListofIntersectionTriangles.h:46
bamg::ListofIntersectionTriangles::Length
double Length()
Definition: ListofIntersectionTriangles.cpp:31
bamg::ListofIntersectionTriangles::SegInterpolation::lBegin
double lBegin
Definition: ListofIntersectionTriangles.h:29
bamg::Norme2_2
RR Norme2_2(const P2< R, RR > x)
Definition: R2.h:89
bamg::AdjacentTriangle::EdgeVertex
BamgVertex * EdgeVertex(const int &i) const
Definition: AdjacentTriangle.cpp:32
bamg::ListofIntersectionTriangles::SplitEdge
void SplitEdge(Mesh &, const R2 &, const R2 &, int nbegin=0)
Definition: ListofIntersectionTriangles.cpp:190
bamg::P2< double, double >
bamg::Mesh::TriangleFindFromCoord
Triangle * TriangleFindFromCoord(const I2 &, long long[3], Triangle *tstart=0)
Definition: Mesh.cpp:3945
bamg::ListofIntersectionTriangles::MaxSize
int MaxSize
Definition: ListofIntersectionTriangles.h:44
bamg::ListofIntersectionTriangles::SegInterpolation::F
R2 F(double s)
Definition: ListofIntersectionTriangles.h:33
bamg::BamgVertex::m
Metric m
Definition: BamgVertex.h:22
bamg::OppositeVertex
static const short OppositeVertex[3]
Definition: macros.h:15
bamg::ListofIntersectionTriangles::NbSeg
int NbSeg
Definition: ListofIntersectionTriangles.h:49
bamgobjects.h
bamg::ListofIntersectionTriangles::SegInterpolation::last
int last
Definition: ListofIntersectionTriangles.h:30
bamg::ListofIntersectionTriangles::ReShape
void ReShape()
Definition: ListofIntersectionTriangles.cpp:175
bamg::NextEdge
static const short NextEdge[3]
Definition: macros.h:17
Generated on Thu Jul 2 2020 08:09:16 for Ice Sheet System Model by   doxygen 1.8.19