Ice Sheet System Model  4.18
Code documentation
Triangle.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  Triangle::Triangle(void){/*{{{*/
14 
15  }
16  /*}}}*/
17  Triangle::Triangle(Mesh *Th,long i,long j,long k) {/*{{{*/
18  BamgVertex *v=Th->vertices;
19  long nbv = Th->nbv;
20  if (i<0 || j<0 || k<0){
21  _error_("i<0 || j<0 || k<0");
22  }
23  if (i>=nbv || j>=nbv || k>=nbv){
24  _error_("i>=nbv || j>=nbv || k>=nbv");
25  }
26  vertices[0]=v+i;
27  vertices[1]=v+j;
28  vertices[2]=v+k;
29  adj[0]=adj[1]=adj[2]=0;
30  AdjEdgeIndex[0]=AdjEdgeIndex[1]=AdjEdgeIndex[2]=0;
31  det=0;
32  }
33  /*}}}*/
34  Triangle::Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2){/*{{{*/
35  vertices[0]=v0;
36  vertices[1]=v1;
37  vertices[2]=v2;
38  adj[0]=adj[1]=adj[2]=0;
39  AdjEdgeIndex[0]=AdjEdgeIndex[1]=AdjEdgeIndex[2]=0;
40  if (v0) det=0;
41  else {
42  det=-1;
43  link=NULL;};
44  }
45  /*}}}*/
46 
47  /*Methods*/
48  AdjacentTriangle Triangle::Adj(int i) const {/*{{{*/
49  return AdjacentTriangle(adj[i],AdjEdgeIndex[i]&3);
50  };/*}}}*/
51  double Triangle::Length() const{/*{{{*/
52 
53  double l;
54 
55  /*Get three vertices A,B and C*/
56  R2 A=*this->vertices[0];
57  R2 B=*this->vertices[1];
58  R2 C=*this->vertices[2];
59 
60  /*Compute edges*/
61  R2 e1=B-A;
62  R2 e2=C-A;
63  R2 e3=B-C;
64 
65  /*Compute edge length*/
66  l=Norme2(e1);
67  l=max(l,Norme2(e2));
68  l=max(l,Norme2(e3));
69 
70  return l;
71  };/*}}}*/
72  void Triangle::Echo(void){/*{{{*/
73 
74  int i;
75 
76  _printf_("Triangle:\n");
77  _printf_(" vertices pointer towards three vertices\n");
78  _printf_(" vertices[0] vertices[1] vertices[2] = " << vertices[0] << " " << vertices[1] << " " << vertices[2] << "\n");
79  _printf_(" adj pointer towards three adjacent triangles\n");
80  _printf_(" adj[0] adj[1] adj[2] = " << adj[0] << " " << adj[1] << " " << adj[2] << "\n");
81  _printf_(" det (integer triangle determinant) = " << det << "\n");
82  if (link){
83  _printf_(" link (pointer toward duplicate triangle)= " << link << "\n");
84  }
85  else{
86  _printf_(" color = " << color << "\n");
87  }
88 
89  _printf_("\nThree vertices:\n");
90  for(i=0;i<3;i++){
91  if (vertices[i]){
92  vertices[i]->Echo();
93  }
94  else{
95  _printf_(" vertex " << i+1 << " does not exist\n");
96  }
97  }
98 
99  return;
100  }
101  /*}}}*/
102  int Triangle::GetAllflag(int a){/*{{{*/
103  return AdjEdgeIndex[a] & 1020;
104  }/*}}}*/
105  int Triangle::Hidden(int a)const {/*{{{*/
106  return AdjEdgeIndex[a]&16;
107  } /*}}}*/
108  int Triangle::Locked(int a)const {/*{{{*/
109  return AdjEdgeIndex[a]&4;
110  } /*}}}*/
111  short Triangle::NuEdgeTriangleAdj(int i) const {/*{{{*/
112  /*Number of the adjacent edge in adj tria (make sure it is between 0 and 2*/
113  return AdjEdgeIndex[i&3]&3;
114  }/*}}}*/
115  long Triangle::Optim(short i,int koption) {/*{{{*/
116  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Optim)*/
117 
118  // turn around (positive direction)
119  Triangle *t=this;
120  long NbSwap =0;
121  int k = 0;
122  int j = OppositeEdge[i];
123  int jp= PreviousEdge[j];
124 
125  // initialize tp, jp the previous triangle & edge
126  Triangle *tp=adj[jp];
127  jp = AdjEdgeIndex[jp]&3;
128  do {
129  while (t->swap(j,koption)){
130  if (k>=20000) _error_("k>=20000");
131  NbSwap++;
132  k++;
133  t= tp->adj[jp]; // set unchange t qnd j for previous triangles
134  j= NextEdge[tp->AdjEdgeIndex[jp]&3];
135  }
136  // end on this Triangle
137  tp = t;
138  jp = NextEdge[j];
139 
140  t= tp->adj[jp]; // set unchange t qnd j for previous triangles
141  j= NextEdge[tp->AdjEdgeIndex[jp]&3];
142 
143  } while( t != this);
144  return NbSwap;
145  }
146  /*}}}*/
147  void Triangle::Renumbering(Triangle *tb,Triangle *te, long *renu){/*{{{*/
148 
149  if (link >=tb && link <te) link = tb + renu[link -tb];
150  if (adj[0] >=tb && adj[0] <te) adj[0] = tb + renu[adj[0]-tb];
151  if (adj[1] >=tb && adj[1] <te) adj[1] = tb + renu[adj[1]-tb];
152  if (adj[2] >=tb && adj[2] <te) adj[2] = tb + renu[adj[2]-tb];
153  }/*}}}*/
154  void Triangle::Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){/*{{{*/
155  if (vertices[0] >=vb && vertices[0] <ve) vertices[0] = vb + renu[vertices[0]-vb];
156  if (vertices[1] >=vb && vertices[1] <ve) vertices[1] = vb + renu[vertices[1]-vb];
157  if (vertices[2] >=vb && vertices[2] <ve) vertices[2] = vb + renu[vertices[2]-vb];
158  }/*}}}*/
159  void Triangle::Set(const Triangle & rec,const Mesh & Th ,Mesh & ThNew){ /*{{{*/
160  *this = rec;
161  if ( vertices[0] ) vertices[0] = ThNew.vertices + Th.GetId(vertices[0]);
162  if ( vertices[1] ) vertices[1] = ThNew.vertices + Th.GetId(vertices[1]);
163  if ( vertices[2] ) vertices[2] = ThNew.vertices + Th.GetId(vertices[2]);
164  if(adj[0]) adj[0] = ThNew.triangles + Th.GetId(adj[0]);
165  if(adj[1]) adj[1] = ThNew.triangles + Th.GetId(adj[1]);
166  if(adj[2]) adj[2] = ThNew.triangles + Th.GetId(adj[2]);
167  if (link >= Th.triangles && link < Th.triangles + Th.nbt)
168  link = ThNew.triangles + Th.GetId(link);
169  }
170  /*}}}*/
171  void Triangle::SetAdjAdj(short a){/*{{{*/
172  // Copy all the mark
173  a &= 3;
174  Triangle *tt=adj[a];
175  AdjEdgeIndex [a] &= 55; // remove MarkUnSwap
176  short aatt = AdjEdgeIndex[a] & 3;
177  if(tt){
178  tt->adj[aatt]=this;
179  tt->AdjEdgeIndex[aatt]=a + (AdjEdgeIndex[a] & 60 ) ;
180  }
181  }/*}}}*/
182  void Triangle::SetAdj2(short a,Triangle *t,short aat){/*{{{*/
183  /*For current triangle:
184  * - a is the index of the edge were the adjency is set (in [0 2])
185  * - t is the adjacent triangle
186  * - aat is the index of the same edge in the adjacent triangle*/
187  adj[a]=t;
188  AdjEdgeIndex[a]=aat;
189  if(t){ //if t!=NULL add adjacent triangle to t (this)
190  t->adj[aat]=this;
191  t->AdjEdgeIndex[aat]=a;
192  }
193  }/*}}}*/
194  void Triangle::SetHidden(int a){/*{{{*/
195  //Get Adjacent Triangle number a
196  Triangle* t = adj[a];
197  //if it exist
198  //C|=D -> C=(C|D) bitwise inclusive OR
199  if(t) t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=16;
200  AdjEdgeIndex[a] |= 16;
201  }/*}}}*/
202  void Triangle::SetLocked(int a){/*{{{*/
203  //mark the edge as on Boundary
204  Triangle * t = adj[a];
205  t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=4;
206  AdjEdgeIndex[a] |= 4;
207  }/*}}}*/
208  void Triangle::SetMarkUnSwap(int a){/*{{{*/
209  Triangle * t = adj[a];
210  t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=8;
211  AdjEdgeIndex[a] |=8 ;
212  }/*}}}*/
213  void Triangle::SetSingleVertexToTriangleConnectivity() { /*{{{*/
214  if (vertices[0]) (vertices[0]->t=this,vertices[0]->IndexInTriangle=0);
215  if (vertices[1]) (vertices[1]->t=this,vertices[1]->IndexInTriangle=1);
216  if (vertices[2]) (vertices[2]->t=this,vertices[2]->IndexInTriangle=2);
217  }/*}}}*/
218  void Triangle::SetUnMarkUnSwap(int a){ /*{{{*/
219  Triangle * t = adj[a];
220  t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] &=55; // 23 + 32
221  AdjEdgeIndex[a] &=55 ;
222  }/*}}}*/
223  Triangle* Triangle::TriangleAdj(int i) const {/*{{{*/
224  return adj[i&3];
225  }/*}}}*/
226  int Triangle::swap(short a,int koption){/*{{{*/
227  /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)*/
228 
229  if(a/4 !=0) return 0;// arete lock or MarkUnSwap
230 
231  Triangle *t1=this,*t2=adj[a];// les 2 triangles adjacent
232  short a1=a,a2=AdjEdgeIndex[a];// les 2 numero de l arete dans les 2 triangles
233  if(a2/4 !=0) return 0; // arete lock or MarkUnSwap
234 
235  BamgVertex *sa=t1->vertices[VerticesOfTriangularEdge[a1][0]];
236  BamgVertex *sb=t1->vertices[VerticesOfTriangularEdge[a1][1]];
237  BamgVertex *s1=t1->vertices[OppositeVertex[a1]];
238  BamgVertex *s2=t2->vertices[OppositeVertex[a2]];
239 
240  long long det1=t1->det , det2=t2->det ;
241  long long detT = det1+det2;
242  long long detA = Abs(det1) + Abs(det2);
243  long long detMin = Min(det1,det2);
244 
245  int OnSwap = 0;
246  // si 2 triangle infini (bord) => detT = -2;
247  if (sa == 0) {// les deux triangles sont frontieres
248  det2=bamg::det(s2->i,sb->i,s1->i);
249  OnSwap = det2 >0;}
250  else if (sb == 0) { // les deux triangles sont frontieres
251  det1=bamg::det(s1->i,sa->i,s2->i);
252  OnSwap = det1 >0;}
253  else if(( s1 != 0) && (s2 != 0) ) {
254  det1 = bamg::det(s1->i,sa->i,s2->i);
255  det2 = detT - det1;
256  OnSwap = (Abs(det1) + Abs(det2)) < detA;
257 
258  long long detMinNew=Min(det1,det2);
259  if (! OnSwap &&(detMinNew>0)) {
260  OnSwap = detMin ==0;
261  if (! OnSwap) {
262  int kopt = koption;
263  while (1)
264  if(kopt) {
265  // critere de Delaunay pure isotrope
266  long long xb1 = sb->i.x - s1->i.x,
267  x21 = s2->i.x - s1->i.x,
268  yb1 = sb->i.y - s1->i.y,
269  y21 = s2->i.y - s1->i.y,
270  xba = sb->i.x - sa->i.x,
271  x2a = s2->i.x - sa->i.x,
272  yba = sb->i.y - sa->i.y,
273  y2a = s2->i.y - sa->i.y;
274  double
275  cosb12 = double(xb1*x21 + yb1*y21),
276  cosba2 = double(xba*x2a + yba*y2a) ,
277  sinb12 = double(det2),
278  sinba2 = double(t2->det);
279 
280  // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2)
281  OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12);
282  break;
283  }
284  else {
285  // critere de Delaunay anisotrope
286  double som;
287  I2 AB=(I2) *sb - (I2) *sa;
288  I2 MAB2=((I2) *sb + (I2) *sa);
289  R2 MAB(MAB2.x*0.5,MAB2.y*0.5);
290  I2 A1=(I2) *s1 - (I2) *sa;
291  I2 D = (I2) * s1 - (I2) * sb ;
292  R2 S2(s2->i.x,s2->i.y);
293  R2 S1(s1->i.x,s1->i.y);
294  {
295  Metric M=s1->m;
296  R2 ABo = M.Orthogonal(AB);
297  R2 A1o = M.Orthogonal(A1);
298  // (A+B)+ x ABo = (S1+B)/2+ y A1
299  // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
300  double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
301  double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
302  if (Abs(d) > dd*1.e-3) {
303  R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
304  som = M.Length(C.x-S2.x,C.y-S2.y) / M.Length(C.x-S1.x,C.y-S1.y);
305  } else
306  {kopt=1;continue;}
307 
308  }
309  {
310  Metric M=s2->m;
311  R2 ABo = M.Orthogonal(AB);
312  R2 A1o = M.Orthogonal(A1);
313  // (A+B)+ x ABo = (S1+B)/2+ y A1
314  // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
315  double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
316  double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
317  if(Abs(d) > dd*1.e-3) {
318  R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
319  som += M.Length(C.x-S2.x,C.y-S2.y) / M.Length(C.x-S1.x,C.y-S1.y);
320  } else
321  {kopt=1;continue;}
322  }
323  OnSwap = som < 2;
324  break;
325  }
326 
327  } // OnSwap
328  } // (! OnSwap &&(det1 > 0) && (det2 > 0) )
329  }
330  if( OnSwap )
331  bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2);
332  else {
333  t1->SetMarkUnSwap(a1);
334  }
335  return OnSwap;
336  }
337  /*}}}*/
338 
339 }
bamg::Mesh::GetId
long GetId(const Triangle &t) const
Definition: Mesh.cpp:2608
bamg::BamgVertex
Definition: BamgVertex.h:15
bamg::BamgVertex::i
I2 i
Definition: BamgVertex.h:20
bamg::Triangle::SetSingleVertexToTriangleConnectivity
void SetSingleVertexToTriangleConnectivity()
Definition: Triangle.cpp:213
bamg::Triangle::Hidden
int Hidden(int a) const
Definition: Triangle.cpp:105
bamg::PreviousEdge
static const short PreviousEdge[3]
Definition: macros.h:18
bamg::Triangle::SetAdjAdj
void SetAdjAdj(short a)
Definition: Triangle.cpp:171
bamg::Triangle::GetAllflag
int GetAllflag(int a)
Definition: Triangle.cpp:102
bamg::Triangle::Optim
long Optim(short a, int=0)
Definition: Triangle.cpp:115
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
bamg
Definition: AdjacentTriangle.cpp:9
bamg::Triangle::SetUnMarkUnSwap
void SetUnMarkUnSwap(int a)
Definition: Triangle.cpp:218
bamg::Triangle
Definition: Triangle.h:13
bamg::Triangle::Locked
int Locked(int a) const
Definition: Triangle.cpp:108
bamg::Metric::Orthogonal
R2 Orthogonal(const R2 x)
Definition: Metric.h:41
bamg::Triangle::link
Triangle * link
Definition: Triangle.h:25
bamg::AdjacentTriangle
Definition: AdjacentTriangle.h:12
bamg::Mesh::triangles
Triangle * triangles
Definition: Mesh.h:28
bamg::OppositeEdge
static const short OppositeEdge[3]
Definition: macros.h:16
bamg::Triangle::Adj
AdjacentTriangle Adj(int i) const
Definition: Triangle.cpp:48
bamg::Triangle::swap
int swap(short a1, int=0)
Definition: Triangle.cpp:226
bamg::Triangle::vertices
BamgVertex * vertices[3]
Definition: Triangle.h:18
bamg::Triangle::SetLocked
void SetLocked(int a)
Definition: Triangle.cpp:202
bamg::Triangle::color
long color
Definition: Triangle.h:26
bamg::Triangle::NuEdgeTriangleAdj
short NuEdgeTriangleAdj(int i) const
Definition: Triangle.cpp:111
det.h
bamg::det
long long det(const I2 &a, const I2 &b, const I2 &c)
Definition: det.h:8
bamg::Norme2
RR Norme2(const P2< R, RR > x)
Definition: R2.h:93
bamg::I2
P2< int, long long > I2
Definition: typedefs.h:11
bamg::Triangle::Length
double Length() const
Definition: Triangle.cpp:51
bamg::Metric
Definition: Metric.h:17
bamg::Mesh::nbt
long nbt
Definition: Mesh.h:35
bamg::Triangle::adj
Triangle * adj[3]
Definition: Triangle.h:19
bamg::BamgVertex::Echo
void Echo()
Definition: BamgVertex.cpp:18
bamg::Mesh
Definition: Mesh.h:21
bamg::Triangle::AdjacentTriangle
friend class AdjacentTriangle
Definition: Triangle.h:15
bamg::Metric::Length
double Length(double Ax, double Ay) const
Definition: Metric.cpp:151
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
bamg::Triangle::det
long long det
Definition: Triangle.h:23
bamg::VerticesOfTriangularEdge
static const short VerticesOfTriangularEdge[3][2]
Definition: macros.h:13
bamg::P2::x
R x
Definition: R2.h:15
bamg::Triangle::AdjEdgeIndex
short AdjEdgeIndex[3]
Definition: Triangle.h:20
bamg::P2::y
R y
Definition: R2.h:15
bamg::Mesh::nbv
long nbv
Definition: Mesh.h:35
bamg::Triangle::SetMarkUnSwap
void SetMarkUnSwap(int a)
Definition: Triangle.cpp:208
bamg::P2< double, double >
bamg::Triangle::TriangleAdj
Triangle * TriangleAdj(int i) const
Definition: Triangle.cpp:223
bamg::Abs
T Abs(const T &a)
Definition: Abs.h:5
bamg::swap
void swap(Triangle *t1, short a1, Triangle *t2, short a2, BamgVertex *s1, BamgVertex *s2, long long det1, long long det2)
Definition: Mesh.cpp:4966
bamg::Triangle::Set
void Set(const Triangle &, const Mesh &, Mesh &)
Definition: Triangle.cpp:159
bamg::Triangle::Renumbering
void Renumbering(Triangle *tb, Triangle *te, long *renu)
Definition: Triangle.cpp:147
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
bamg::Min
T Min(const T &a, const T &b)
Definition: extrema.h:6
bamg::Mesh::vertices
BamgVertex * vertices
Definition: Mesh.h:27
bamg::Triangle::Triangle
Triangle()
Definition: Triangle.cpp:13
bamg::BamgVertex::m
Metric m
Definition: BamgVertex.h:22
bamg::OppositeVertex
static const short OppositeVertex[3]
Definition: macros.h:15
bamg::Triangle::Echo
void Echo()
Definition: Triangle.cpp:72
bamg::Triangle::SetAdj2
void SetAdj2(short a, Triangle *t, short aat)
Definition: Triangle.cpp:182
bamgobjects.h
bamg::Triangle::SetHidden
void SetHidden(int a)
Definition: Triangle.cpp:194
bamg::BamgVertex::IndexInTriangle
short IndexInTriangle
Definition: BamgVertex.h:25
bamg::NextEdge
static const short NextEdge[3]
Definition: macros.h:17
bamg::BamgVertex::t
Triangle * t
Definition: BamgVertex.h:28
Generated on Thu Jul 2 2020 08:09:16 for Ice Sheet System Model by   doxygen 1.8.19