[3913] | 1 | #include <cstdio>
|
---|
| 2 | #include <cstring>
|
---|
| 3 | #include <cmath>
|
---|
| 4 | #include <ctime>
|
---|
| 5 |
|
---|
[13530] | 6 | #include "../objects/objects.h"
|
---|
[14949] | 7 | #include "./det.h"
|
---|
[3913] | 8 |
|
---|
| 9 | namespace bamg {
|
---|
| 10 |
|
---|
| 11 | /*Constructors/Destructors*/
|
---|
[12365] | 12 | /*FUNCTION Triangle(){{{*/
|
---|
[9371] | 13 | Triangle::Triangle(void){
|
---|
| 14 |
|
---|
| 15 | }
|
---|
| 16 | /*}}}*/
|
---|
[12365] | 17 | /*FUNCTION Triangle(Mesh *Th,long i,long j,long k) {{{*/
|
---|
[5095] | 18 | Triangle::Triangle(Mesh *Th,long i,long j,long k) {
|
---|
[5120] | 19 | BamgVertex *v=Th->vertices;
|
---|
[3913] | 20 | long nbv = Th->nbv;
|
---|
| 21 | if (i<0 || j<0 || k<0){
|
---|
[13036] | 22 | _error_("i<0 || j<0 || k<0");
|
---|
[3913] | 23 | }
|
---|
| 24 | if (i>=nbv || j>=nbv || k>=nbv){
|
---|
[13036] | 25 | _error_("i>=nbv || j>=nbv || k>=nbv");
|
---|
[3913] | 26 | }
|
---|
[5148] | 27 | vertices[0]=v+i;
|
---|
| 28 | vertices[1]=v+j;
|
---|
| 29 | vertices[2]=v+k;
|
---|
| 30 | adj[0]=adj[1]=adj[2]=0;
|
---|
[5401] | 31 | AdjEdgeIndex[0]=AdjEdgeIndex[1]=AdjEdgeIndex[2]=0;
|
---|
[3913] | 32 | det=0;
|
---|
| 33 | }
|
---|
| 34 | /*}}}*/
|
---|
[12365] | 35 | /*FUNCTION Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2) {{{*/
|
---|
[5120] | 36 | Triangle::Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2){
|
---|
[5148] | 37 | vertices[0]=v0;
|
---|
| 38 | vertices[1]=v1;
|
---|
| 39 | vertices[2]=v2;
|
---|
| 40 | adj[0]=adj[1]=adj[2]=0;
|
---|
[5401] | 41 | AdjEdgeIndex[0]=AdjEdgeIndex[1]=AdjEdgeIndex[2]=0;
|
---|
[3913] | 42 | if (v0) det=0;
|
---|
| 43 | else {
|
---|
| 44 | det=-1;
|
---|
| 45 | link=NULL;};
|
---|
| 46 | }
|
---|
| 47 | /*}}}*/
|
---|
| 48 |
|
---|
| 49 | /*Methods*/
|
---|
[12365] | 50 | /*FUNCTION Triangle::Adj{{{*/
|
---|
[5151] | 51 | AdjacentTriangle Triangle::Adj(int i) const {
|
---|
[5401] | 52 | return AdjacentTriangle(adj[i],AdjEdgeIndex[i]&3);
|
---|
[5151] | 53 | };/*}}}*/
|
---|
[12365] | 54 | /*FUNCTION Triangle::Anisotropy{{{*/
|
---|
[11584] | 55 | double Triangle::Anisotropy() const{
|
---|
| 56 |
|
---|
| 57 | double lmin,lmax;
|
---|
| 58 |
|
---|
| 59 | /*Get three vertices A,B and C*/
|
---|
| 60 | R2 A=*this->vertices[0];
|
---|
| 61 | R2 B=*this->vertices[1];
|
---|
| 62 | R2 C=*this->vertices[2];
|
---|
| 63 |
|
---|
| 64 | /*Compute edges*/
|
---|
| 65 | R2 e1=B-A;
|
---|
| 66 | R2 e2=C-A;
|
---|
| 67 | R2 e3=B-C;
|
---|
| 68 |
|
---|
| 69 | /*Compute edge length*/
|
---|
| 70 | double l1=Norme2(e1);
|
---|
| 71 | double l2=Norme2(e2);
|
---|
| 72 | double l3=Norme2(e3);
|
---|
| 73 |
|
---|
| 74 | lmin=l1;
|
---|
| 75 | lmin=min(lmin,l2);
|
---|
| 76 | lmin=min(lmin,l3);
|
---|
| 77 | lmax=l1;
|
---|
| 78 | lmax=max(lmax,l2);
|
---|
| 79 | lmax=max(lmax,l3);
|
---|
| 80 |
|
---|
| 81 | return lmax/lmin;
|
---|
| 82 | };/*}}}*/
|
---|
[12365] | 83 | /*FUNCTION Triangle::Length{{{*/
|
---|
[11584] | 84 | double Triangle::Length() const{
|
---|
| 85 |
|
---|
| 86 | double l;
|
---|
| 87 |
|
---|
| 88 | /*Get three vertices A,B and C*/
|
---|
| 89 | R2 A=*this->vertices[0];
|
---|
| 90 | R2 B=*this->vertices[1];
|
---|
| 91 | R2 C=*this->vertices[2];
|
---|
| 92 |
|
---|
| 93 | /*Compute edges*/
|
---|
| 94 | R2 e1=B-A;
|
---|
| 95 | R2 e2=C-A;
|
---|
| 96 | R2 e3=B-C;
|
---|
| 97 |
|
---|
| 98 | /*Compute edge length*/
|
---|
| 99 | l=Norme2(e1);
|
---|
| 100 | l=max(l,Norme2(e2));
|
---|
| 101 | l=max(l,Norme2(e3));
|
---|
| 102 |
|
---|
| 103 | return l;
|
---|
| 104 | };/*}}}*/
|
---|
[12365] | 105 | /*FUNCTION Triangle::Echo {{{*/
|
---|
[3913] | 106 | void Triangle::Echo(void){
|
---|
| 107 |
|
---|
| 108 | int i;
|
---|
| 109 |
|
---|
[12509] | 110 | _printLine_("Triangle:");
|
---|
| 111 | _printLine_(" vertices pointer towards three vertices");
|
---|
| 112 | _printLine_(" vertices[0] vertices[1] vertices[2] = " << vertices[0] << " " << vertices[1] << " " << vertices[2]);
|
---|
| 113 | _printLine_(" adj pointer towards three adjacent triangles");
|
---|
| 114 | _printLine_(" adj[0] adj[1] adj[2] = " << adj[0] << " " << adj[1] << " " << adj[2]);
|
---|
| 115 | _printLine_(" det (integer triangle determinant) = " << det);
|
---|
[3913] | 116 | if (link){
|
---|
[12509] | 117 | _printLine_(" link (pointer toward duplicate triangle)= " << link);
|
---|
[3913] | 118 | }
|
---|
| 119 | else{
|
---|
[12509] | 120 | _printLine_(" color = " << color);
|
---|
[3913] | 121 | }
|
---|
| 122 |
|
---|
[12509] | 123 | _printLine_("\nThree vertices:");
|
---|
[3913] | 124 | for(i=0;i<3;i++){
|
---|
[5148] | 125 | if (vertices[i]){
|
---|
| 126 | vertices[i]->Echo();
|
---|
[3913] | 127 | }
|
---|
| 128 | else{
|
---|
[12509] | 129 | _printLine_(" vertex " << i+1 << " does not exist");
|
---|
[3913] | 130 | }
|
---|
| 131 | }
|
---|
| 132 |
|
---|
| 133 | return;
|
---|
| 134 | }
|
---|
| 135 | /*}}}*/
|
---|
[12365] | 136 | /*FUNCTION Triangle::FindBoundaryEdge{{{*/
|
---|
[5143] | 137 | AdjacentTriangle Triangle::FindBoundaryEdge(int i) const{
|
---|
[3913] | 138 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindBoundaryEdge)*/
|
---|
| 139 |
|
---|
| 140 | /*Intermediary*/
|
---|
| 141 | Triangle* ttc=NULL;
|
---|
| 142 | int k,j,jc;
|
---|
| 143 |
|
---|
| 144 | // call current triangle t
|
---|
| 145 | Triangle* t = (Triangle*)this;
|
---|
| 146 |
|
---|
| 147 | //is the current triangle inside or outside?
|
---|
| 148 | int outside=!link ;
|
---|
| 149 |
|
---|
| 150 | // EdgesVertexTriangle[3][2] = {{1,2},{2,0},{0,1}};
|
---|
| 151 | // initialize j as the first vertex of the ith edge
|
---|
| 152 | j=EdgesVertexTriangle[i][0];
|
---|
| 153 |
|
---|
| 154 | //Loop over the adjacent triangle of t
|
---|
| 155 | k=0;
|
---|
| 156 | do{
|
---|
| 157 | //keep track of outside
|
---|
| 158 | int outsidep = outside;
|
---|
| 159 | //increment k
|
---|
| 160 | k++;
|
---|
| 161 | //Get ttc, adjacent triangle of t with respect to vertex j
|
---|
[5148] | 162 | ttc = t->adj[j];
|
---|
[3913] | 163 | //is the current triangle inside or outside?
|
---|
| 164 | outside = !ttc->link;
|
---|
| 165 | //if both previous triangle are outside, return
|
---|
[5143] | 166 | if (outside+outsidep == 1) return AdjacentTriangle(t,j);
|
---|
[3913] | 167 |
|
---|
| 168 | //update t and j
|
---|
| 169 | t = ttc;
|
---|
| 170 | //NextEdge[3] = {1,2,0};
|
---|
[5401] | 171 | jc = NextEdge[t->AdjEdgeIndex[j]&3];
|
---|
[3913] | 172 | j = NextEdge[jc];
|
---|
| 173 |
|
---|
| 174 | //check number of iterations
|
---|
| 175 | if (k>=2000){
|
---|
[13036] | 176 | _error_("too many iteration in Triangle::FindBoundaryEdge (k>=2000)");
|
---|
[3913] | 177 | }
|
---|
| 178 | } while (this!= t);
|
---|
| 179 | //not found, return empty triangle
|
---|
[5143] | 180 | return AdjacentTriangle(NULL,0);
|
---|
[3913] | 181 | }
|
---|
[12365] | 182 | /*}}}*/
|
---|
| 183 | /*FUNCTION Triangle::GetAllflag{{{*/
|
---|
[5151] | 184 | int Triangle::GetAllflag(int a){
|
---|
[5401] | 185 | return AdjEdgeIndex[a] & 1020;
|
---|
[5151] | 186 | }/*}}}*/
|
---|
[12365] | 187 | /*FUNCTION Triangle::Hidden{{{*/
|
---|
[5151] | 188 | int Triangle::Hidden(int a)const {
|
---|
[5401] | 189 | return AdjEdgeIndex[a]&16;
|
---|
[5151] | 190 | } /*}}}*/
|
---|
[12365] | 191 | /*FUNCTION Triangle::Locked{{{*/
|
---|
[5151] | 192 | int Triangle::Locked(int a)const {
|
---|
[5401] | 193 | return AdjEdgeIndex[a]&4;
|
---|
[5151] | 194 | } /*}}}*/
|
---|
[12365] | 195 | /*FUNCTION Triangle::NuEdgeTriangleAdj{{{*/
|
---|
[5151] | 196 | short Triangle::NuEdgeTriangleAdj(int i) const {
|
---|
[5187] | 197 | /*Number of the adjacent edge in adj tria (make sure it is between 0 and 2*/
|
---|
[5401] | 198 | return AdjEdgeIndex[i&3]&3;
|
---|
[5151] | 199 | }/*}}}*/
|
---|
[12365] | 200 | /*FUNCTION Triangle::Optim{{{*/
|
---|
[3913] | 201 | long Triangle::Optim(short i,int koption) {
|
---|
| 202 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Optim)*/
|
---|
| 203 |
|
---|
| 204 | // turn around (positive direction)
|
---|
| 205 | Triangle *t=this;
|
---|
| 206 | long NbSwap =0;
|
---|
| 207 | int k = 0;
|
---|
| 208 | int j = OppositeEdge[i];
|
---|
| 209 | int jp= PreviousEdge[j];
|
---|
| 210 |
|
---|
| 211 | // initialize tp, jp the previous triangle & edge
|
---|
[5148] | 212 | Triangle *tp=adj[jp];
|
---|
[5401] | 213 | jp = AdjEdgeIndex[jp]&3;
|
---|
[3913] | 214 | do {
|
---|
| 215 | while (t->swap(j,koption)){
|
---|
[13036] | 216 | if (k>=20000) _error_("k>=20000");
|
---|
[3913] | 217 | NbSwap++;
|
---|
| 218 | k++;
|
---|
[5148] | 219 | t= tp->adj[jp]; // set unchange t qnd j for previous triangles
|
---|
[5401] | 220 | j= NextEdge[tp->AdjEdgeIndex[jp]&3];
|
---|
[3913] | 221 | }
|
---|
| 222 | // end on this Triangle
|
---|
| 223 | tp = t;
|
---|
| 224 | jp = NextEdge[j];
|
---|
| 225 |
|
---|
[5148] | 226 | t= tp->adj[jp]; // set unchange t qnd j for previous triangles
|
---|
[5401] | 227 | j= NextEdge[tp->AdjEdgeIndex[jp]&3];
|
---|
[3913] | 228 |
|
---|
| 229 | } while( t != this);
|
---|
| 230 | return NbSwap;
|
---|
| 231 | }
|
---|
[12365] | 232 | /*}}}*/
|
---|
| 233 | /*FUNCTION Triangle::Quadrangle {{{*/
|
---|
[5120] | 234 | Triangle* Triangle::Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const{
|
---|
[3913] | 235 | // return the other triangle of the quad if a quad or 0 if not a quat
|
---|
| 236 | Triangle * t =0;
|
---|
| 237 | if (link) {
|
---|
| 238 | int a=-1;
|
---|
[5401] | 239 | if (AdjEdgeIndex[0] & 16 ) a=0;
|
---|
| 240 | if (AdjEdgeIndex[1] & 16 ) a=1;
|
---|
| 241 | if (AdjEdgeIndex[2] & 16 ) a=2;
|
---|
[3913] | 242 | if (a>=0) {
|
---|
[5148] | 243 | t = adj[a];
|
---|
[3913] | 244 | // if (t-this<0) return 0;
|
---|
[5148] | 245 | v2 = vertices[VerticesOfTriangularEdge[a][0]];
|
---|
| 246 | v0 = vertices[VerticesOfTriangularEdge[a][1]];
|
---|
| 247 | v1 = vertices[OppositeEdge[a]];
|
---|
[5401] | 248 | v3 = t->vertices[OppositeEdge[AdjEdgeIndex[a]&3]];
|
---|
[3913] | 249 | }
|
---|
| 250 | }
|
---|
| 251 | return t;
|
---|
| 252 | }
|
---|
| 253 | /*}}}*/
|
---|
[12365] | 254 | /*FUNCTION Triangle::QualityQuad {{{*/
|
---|
[3913] | 255 | double Triangle::QualityQuad(int a,int option) const{
|
---|
| 256 | double q;
|
---|
[5401] | 257 | if (!link || AdjEdgeIndex[a] &4)
|
---|
[3913] | 258 | q= -1;
|
---|
| 259 | else {
|
---|
[5148] | 260 | Triangle * t = adj[a];
|
---|
[3913] | 261 | if (t-this<0) q= -1;// because we do 2 times
|
---|
| 262 | else if (!t->link ) q= -1;
|
---|
[5401] | 263 | else if (AdjEdgeIndex[0] & 16 || AdjEdgeIndex[1] & 16 || AdjEdgeIndex[2] & 16 || t->AdjEdgeIndex[0] & 16 || t->AdjEdgeIndex[1] & 16 || t->AdjEdgeIndex[2] & 16 )
|
---|
[3913] | 264 | q= -1;
|
---|
| 265 | else if(option){
|
---|
[5148] | 266 | const BamgVertex & v2 = *vertices[VerticesOfTriangularEdge[a][0]];
|
---|
| 267 | const BamgVertex & v0 = *vertices[VerticesOfTriangularEdge[a][1]];
|
---|
| 268 | const BamgVertex & v1 = *vertices[OppositeEdge[a]];
|
---|
[5401] | 269 | const BamgVertex & v3 = * t->vertices[OppositeEdge[AdjEdgeIndex[a]&3]];
|
---|
[3913] | 270 | q = QuadQuality(v0,v1,v2,v3); // do the float part
|
---|
| 271 | }
|
---|
| 272 | else q= 1;
|
---|
| 273 | }
|
---|
| 274 | return q;
|
---|
| 275 | }
|
---|
| 276 | /*}}}*/
|
---|
[12365] | 277 | /*FUNCTION Triangle::Renumbering(Triangle *tb,Triangle *te, long *renu){{{*/
|
---|
[5151] | 278 | void Triangle::Renumbering(Triangle *tb,Triangle *te, long *renu){
|
---|
| 279 |
|
---|
| 280 | if (link >=tb && link <te) link = tb + renu[link -tb];
|
---|
| 281 | if (adj[0] >=tb && adj[0] <te) adj[0] = tb + renu[adj[0]-tb];
|
---|
| 282 | if (adj[1] >=tb && adj[1] <te) adj[1] = tb + renu[adj[1]-tb];
|
---|
| 283 | if (adj[2] >=tb && adj[2] <te) adj[2] = tb + renu[adj[2]-tb];
|
---|
| 284 | }/*}}}*/
|
---|
[12365] | 285 | /*FUNCTION Triangle::Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){{{*/
|
---|
[5151] | 286 | void Triangle::Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
|
---|
| 287 | if (vertices[0] >=vb && vertices[0] <ve) vertices[0] = vb + renu[vertices[0]-vb];
|
---|
| 288 | if (vertices[1] >=vb && vertices[1] <ve) vertices[1] = vb + renu[vertices[1]-vb];
|
---|
| 289 | if (vertices[2] >=vb && vertices[2] <ve) vertices[2] = vb + renu[vertices[2]-vb];
|
---|
| 290 | }/*}}}*/
|
---|
[12365] | 291 | /*FUNCTION Triangle::Set {{{*/
|
---|
[5095] | 292 | void Triangle::Set(const Triangle & rec,const Mesh & Th ,Mesh & ThNew){
|
---|
[3913] | 293 | *this = rec;
|
---|
[5149] | 294 | if ( vertices[0] ) vertices[0] = ThNew.vertices + Th.GetId(vertices[0]);
|
---|
| 295 | if ( vertices[1] ) vertices[1] = ThNew.vertices + Th.GetId(vertices[1]);
|
---|
| 296 | if ( vertices[2] ) vertices[2] = ThNew.vertices + Th.GetId(vertices[2]);
|
---|
| 297 | if(adj[0]) adj[0] = ThNew.triangles + Th.GetId(adj[0]);
|
---|
| 298 | if(adj[1]) adj[1] = ThNew.triangles + Th.GetId(adj[1]);
|
---|
| 299 | if(adj[2]) adj[2] = ThNew.triangles + Th.GetId(adj[2]);
|
---|
[3913] | 300 | if (link >= Th.triangles && link < Th.triangles + Th.nbt)
|
---|
[5149] | 301 | link = ThNew.triangles + Th.GetId(link);
|
---|
[3913] | 302 | }
|
---|
| 303 | /*}}}*/
|
---|
[12365] | 304 | /*FUNCTION Triangle::SetAdjAdj{{{*/
|
---|
[5151] | 305 | void Triangle::SetAdjAdj(short a){
|
---|
| 306 | // Copy all the mark
|
---|
| 307 | a &= 3;
|
---|
| 308 | register Triangle *tt=adj[a];
|
---|
[5401] | 309 | AdjEdgeIndex [a] &= 55; // remove MarkUnSwap
|
---|
| 310 | register short aatt = AdjEdgeIndex[a] & 3;
|
---|
[5151] | 311 | if(tt){
|
---|
| 312 | tt->adj[aatt]=this;
|
---|
[5401] | 313 | tt->AdjEdgeIndex[aatt]=a + (AdjEdgeIndex[a] & 60 ) ;
|
---|
[5151] | 314 | }
|
---|
| 315 | }/*}}}*/
|
---|
[12365] | 316 | /*FUNCTION Triangle::SetAdj2{{{*/
|
---|
[5151] | 317 | void Triangle::SetAdj2(short a,Triangle *t,short aat){
|
---|
[5605] | 318 | /*For current triangle:
|
---|
| 319 | * - a is the index of the edge were the adjency is set (in [0 2])
|
---|
| 320 | * - t is the adjacent triangle
|
---|
| 321 | * - aat is the index of the same edge in the adjacent triangle*/
|
---|
| 322 | adj[a]=t;
|
---|
| 323 | AdjEdgeIndex[a]=aat;
|
---|
| 324 | if(t){ //if t!=NULL add adjacent triangle to t (this)
|
---|
[5151] | 325 | t->adj[aat]=this;
|
---|
[5401] | 326 | t->AdjEdgeIndex[aat]=a;
|
---|
[5151] | 327 | }
|
---|
| 328 | }/*}}}*/
|
---|
[12365] | 329 | /*FUNCTION Triangle::SetAllFlag{{{*/
|
---|
[5151] | 330 | void Triangle::SetAllFlag(int a,int f){
|
---|
[5401] | 331 | AdjEdgeIndex[a] = (AdjEdgeIndex[a] &3) + (1020 & f);
|
---|
[5151] | 332 | }/*}}}*/
|
---|
[12365] | 333 | /*FUNCTION Triangle::SetDet{{{*/
|
---|
[5151] | 334 | void Triangle::SetDet() {
|
---|
| 335 | if(vertices[0] && vertices[1] && vertices[2]) det = bamg::det(*vertices[0],*vertices[1],*vertices[2]);
|
---|
| 336 | else det = -1;
|
---|
| 337 | }/*}}}*/
|
---|
[12365] | 338 | /*FUNCTION Triangle::SetHidden{{{*/
|
---|
[5151] | 339 | void Triangle::SetHidden(int a){
|
---|
| 340 | //Get Adjacent Triangle number a
|
---|
| 341 | register Triangle* t = adj[a];
|
---|
| 342 | //if it exist
|
---|
| 343 | //C|=D -> C=(C|D) bitwise inclusive OR
|
---|
[5401] | 344 | if(t) t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=16;
|
---|
| 345 | AdjEdgeIndex[a] |= 16;
|
---|
[5151] | 346 | }/*}}}*/
|
---|
[12365] | 347 | /*FUNCTION Triangle::SetLocked{{{*/
|
---|
[5151] | 348 | void Triangle::SetLocked(int a){
|
---|
| 349 | //mark the edge as on Boundary
|
---|
| 350 | register Triangle * t = adj[a];
|
---|
[5401] | 351 | t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=4;
|
---|
| 352 | AdjEdgeIndex[a] |= 4;
|
---|
[5151] | 353 | }/*}}}*/
|
---|
[12365] | 354 | /*FUNCTION Triangle::SetMarkUnSwap{{{*/
|
---|
[5151] | 355 | void Triangle::SetMarkUnSwap(int a){
|
---|
| 356 | register Triangle * t = adj[a];
|
---|
[5401] | 357 | t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] |=8;
|
---|
| 358 | AdjEdgeIndex[a] |=8 ;
|
---|
[5151] | 359 | }/*}}}*/
|
---|
[12365] | 360 | /*FUNCTION Triangle::SetSingleVertexToTriangleConnectivity{{{*/
|
---|
[5151] | 361 | void Triangle::SetSingleVertexToTriangleConnectivity() {
|
---|
[5460] | 362 | if (vertices[0]) (vertices[0]->t=this,vertices[0]->IndexInTriangle=0);
|
---|
| 363 | if (vertices[1]) (vertices[1]->t=this,vertices[1]->IndexInTriangle=1);
|
---|
| 364 | if (vertices[2]) (vertices[2]->t=this,vertices[2]->IndexInTriangle=2);
|
---|
[5151] | 365 | }/*}}}*/
|
---|
[12365] | 366 | /*FUNCTION Triangle::SetUnMarkUnSwap{{{*/
|
---|
[5151] | 367 | void Triangle::SetUnMarkUnSwap(int a){
|
---|
| 368 | register Triangle * t = adj[a];
|
---|
[5401] | 369 | t->AdjEdgeIndex[AdjEdgeIndex[a] & 3] &=55; // 23 + 32
|
---|
| 370 | AdjEdgeIndex[a] &=55 ;
|
---|
[5151] | 371 | }/*}}}*/
|
---|
[12365] | 372 | /*FUNCTION Triangle::swap{{{*/
|
---|
[3913] | 373 | int Triangle::swap(short a,int koption){
|
---|
| 374 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)*/
|
---|
| 375 |
|
---|
| 376 | if(a/4 !=0) return 0;// arete lock or MarkUnSwap
|
---|
| 377 |
|
---|
[5148] | 378 | register Triangle *t1=this,*t2=adj[a];// les 2 triangles adjacent
|
---|
[5401] | 379 | register short a1=a,a2=AdjEdgeIndex[a];// les 2 numero de l arete dans les 2 triangles
|
---|
[3913] | 380 | if(a2/4 !=0) return 0; // arete lock or MarkUnSwap
|
---|
| 381 |
|
---|
[5148] | 382 | register BamgVertex *sa=t1->vertices[VerticesOfTriangularEdge[a1][0]];
|
---|
| 383 | register BamgVertex *sb=t1->vertices[VerticesOfTriangularEdge[a1][1]];
|
---|
| 384 | register BamgVertex *s1=t1->vertices[OppositeVertex[a1]];
|
---|
| 385 | register BamgVertex *s2=t2->vertices[OppositeVertex[a2]];
|
---|
[3913] | 386 |
|
---|
| 387 | Icoor2 det1=t1->det , det2=t2->det ;
|
---|
| 388 | Icoor2 detT = det1+det2;
|
---|
| 389 | Icoor2 detA = Abs(det1) + Abs(det2);
|
---|
| 390 | Icoor2 detMin = Min(det1,det2);
|
---|
| 391 |
|
---|
| 392 | int OnSwap = 0;
|
---|
| 393 | // si 2 triangle infini (bord) => detT = -2;
|
---|
| 394 | if (sa == 0) {// les deux triangles sont frontieres
|
---|
| 395 | det2=bamg::det(s2->i,sb->i,s1->i);
|
---|
| 396 | OnSwap = det2 >0;}
|
---|
| 397 | else if (sb == 0) { // les deux triangles sont frontieres
|
---|
| 398 | det1=bamg::det(s1->i,sa->i,s2->i);
|
---|
| 399 | OnSwap = det1 >0;}
|
---|
| 400 | else if(( s1 != 0) && (s2 != 0) ) {
|
---|
| 401 | det1 = bamg::det(s1->i,sa->i,s2->i);
|
---|
| 402 | det2 = detT - det1;
|
---|
| 403 | OnSwap = (Abs(det1) + Abs(det2)) < detA;
|
---|
| 404 |
|
---|
| 405 | Icoor2 detMinNew=Min(det1,det2);
|
---|
| 406 | // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test
|
---|
| 407 | if (! OnSwap &&(detMinNew>0)) {
|
---|
| 408 | OnSwap = detMin ==0;
|
---|
| 409 | if (! OnSwap) {
|
---|
| 410 | int kopt = koption;
|
---|
| 411 | while (1)
|
---|
| 412 | if(kopt) {
|
---|
| 413 | // critere de Delaunay pure isotrope
|
---|
| 414 | register Icoor2 xb1 = sb->i.x - s1->i.x,
|
---|
| 415 | x21 = s2->i.x - s1->i.x,
|
---|
| 416 | yb1 = sb->i.y - s1->i.y,
|
---|
| 417 | y21 = s2->i.y - s1->i.y,
|
---|
| 418 | xba = sb->i.x - sa->i.x,
|
---|
| 419 | x2a = s2->i.x - sa->i.x,
|
---|
| 420 | yba = sb->i.y - sa->i.y,
|
---|
| 421 | y2a = s2->i.y - sa->i.y;
|
---|
| 422 | register double
|
---|
| 423 | cosb12 = double(xb1*x21 + yb1*y21),
|
---|
| 424 | cosba2 = double(xba*x2a + yba*y2a) ,
|
---|
| 425 | sinb12 = double(det2),
|
---|
| 426 | sinba2 = double(t2->det);
|
---|
| 427 |
|
---|
| 428 | // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2)
|
---|
| 429 | OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12);
|
---|
| 430 | break;
|
---|
| 431 | }
|
---|
| 432 | else {
|
---|
| 433 | // critere de Delaunay anisotrope
|
---|
| 434 | double som;
|
---|
| 435 | I2 AB=(I2) *sb - (I2) *sa;
|
---|
| 436 | I2 MAB2=((I2) *sb + (I2) *sa);
|
---|
| 437 | R2 MAB(MAB2.x*0.5,MAB2.y*0.5);
|
---|
| 438 | I2 A1=(I2) *s1 - (I2) *sa;
|
---|
| 439 | I2 D = (I2) * s1 - (I2) * sb ;
|
---|
| 440 | R2 S2(s2->i.x,s2->i.y);
|
---|
| 441 | R2 S1(s1->i.x,s1->i.y);
|
---|
| 442 | {
|
---|
| 443 | Metric M=s1->m;
|
---|
| 444 | R2 ABo = M.Orthogonal(AB);
|
---|
| 445 | R2 A1o = M.Orthogonal(A1);
|
---|
| 446 | // (A+B)+ x ABo = (S1+B)/2+ y A1
|
---|
| 447 | // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
|
---|
| 448 | double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
|
---|
| 449 | double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
|
---|
| 450 | if (Abs(d) > dd*1.e-3) {
|
---|
| 451 | R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
|
---|
| 452 | som = M(C - S2)/M(C - S1);
|
---|
| 453 | } else
|
---|
| 454 | {kopt=1;continue;}
|
---|
| 455 |
|
---|
| 456 | }
|
---|
| 457 | {
|
---|
| 458 | Metric M=s2->m;
|
---|
| 459 | R2 ABo = M.Orthogonal(AB);
|
---|
| 460 | R2 A1o = M.Orthogonal(A1);
|
---|
| 461 | // (A+B)+ x ABo = (S1+B)/2+ y A1
|
---|
| 462 | // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2
|
---|
| 463 | double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x);
|
---|
| 464 | double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2
|
---|
| 465 | if(Abs(d) > dd*1.e-3) {
|
---|
| 466 | R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d));
|
---|
| 467 | som += M(C - S2)/M(C - S1);
|
---|
| 468 | } else
|
---|
| 469 | {kopt=1;continue;}
|
---|
| 470 | }
|
---|
| 471 | OnSwap = som < 2;
|
---|
| 472 | break;
|
---|
| 473 | }
|
---|
| 474 |
|
---|
| 475 | } // OnSwap
|
---|
| 476 | } // (! OnSwap &&(det1 > 0) && (det2 > 0) )
|
---|
| 477 | }
|
---|
| 478 | if( OnSwap )
|
---|
| 479 | bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2);
|
---|
| 480 | else {
|
---|
| 481 | t1->SetMarkUnSwap(a1);
|
---|
| 482 | }
|
---|
| 483 | return OnSwap;
|
---|
| 484 | }
|
---|
[12365] | 485 | /*}}}*/
|
---|
| 486 | /*FUNCTION Triangle::TriangleAdj{{{*/
|
---|
[5151] | 487 | Triangle* Triangle::TriangleAdj(int i) const {
|
---|
| 488 | return adj[i&3];
|
---|
| 489 | }/*}}}*/
|
---|
[3913] | 490 |
|
---|
| 491 | }
|
---|