![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3239
- Timestamp:
- 03/10/10 07:47:43 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 12 deleted
- 3 edited
-
BamgObjects.h (modified) (1 diff)
-
objects/Triangles.cpp (modified) (1 diff)
-
shared/AGoodNumberPrimeWith.cpp (deleted)
-
shared/AGoodNumberPrimeWith.h (deleted)
-
shared/CloseBoundaryEdge.cpp (deleted)
-
shared/CloseBoundaryEdge.h (deleted)
-
shared/CloseBoundaryEdgeV2.cpp (deleted)
-
shared/CloseBoundaryEdgeV2.h (deleted)
-
shared/ForceEdge.cpp (deleted)
-
shared/ForceEdge.h (deleted)
-
shared/SwapForForcingEdge.cpp (deleted)
-
shared/SwapForForcingEdge.h (deleted)
-
shared/shared.h (modified) (1 diff)
-
shared/swap.cpp (deleted)
-
shared/swap.h (deleted)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/BamgObjects.h
r3238 r3239 334 334 /*}}}1*/ 335 335 336 /*Other prototypes {{{1*/ 336 /*Other prototypes IN TRIANGLES.CPP (TO BE REORGANIZED){{{1*/ 337 Int4 AGoodNumberPrimeWith(Int4 n); 337 338 TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ; 338 339 TriangleAdjacent CloseBoundaryEdgeV2(I2 A,Triangle *t, double &a,double &b); 339 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside);340 340 void swap(Triangle *t1,Int1 a1, 341 341 Triangle *t2,Int1 a2, -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r3238 r3239 5359 5359 /*}}}1*/ 5360 5360 5361 /*Intermediary*/ 5362 /*FUNCTION swap{{{1*/ 5363 void swap(Triangle *t1,Int1 a1, Triangle *t2,Int1 a2, Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2){ 5364 // -------------------------------------------------------------- 5365 // Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles 5366 // 5367 // sb sb 5368 // / | \ / \ ! 5369 // as1/ | \ /a2 \ ! 5370 // / | \ / t2 \ ! 5371 // s1 /t1 | t2 \s2 --> s1 /___as2___\s2 ! 5372 // \ a1|a2 / \ as1 / 5373 // \ | / \ t1 / 5374 // \ | / as2 \ a1/ 5375 // \ | / \ / 5376 // sa sa 5377 // ------------------------------------------------------------- 5378 int as1 = NextEdge[a1]; 5379 int as2 = NextEdge[a2]; 5380 int ap1 = PreviousEdge[a1]; 5381 int ap2 = PreviousEdge[a2]; 5382 (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb 5383 (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa 5384 // mise a jour des 2 adjacences externes 5385 TriangleAdjacent taas1 = t1->Adj(as1), 5386 taas2 = t2->Adj(as2), 5387 tas1(t1,as1), tas2(t2,as2), 5388 ta1(t1,a1),ta2(t2,a2); 5389 // externe haut gauche 5390 taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap()); 5391 // externe bas droite 5392 taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap()); 5393 // remove the Mark UnMarkSwap 5394 t1->SetUnMarkUnSwap(ap1); 5395 t2->SetUnMarkUnSwap(ap2); 5396 // interne 5397 tas1.SetAdj2(tas2); 5398 5399 t1->det = det1; 5400 t2->det = det2; 5401 5402 t1->SetTriangleContainingTheVertex(); 5403 t2->SetTriangleContainingTheVertex(); 5404 } // end swap 5405 /*}}}1*/ 5406 /*FUNCTION FindTriangle{{{1*/ 5407 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside){ 5408 I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); 5409 Icoor2 dete[3]; 5410 Triangle & tb = *Th.FindTriangleContaining(I,dete); 5411 5412 if (tb.link) 5413 { // internal point in a true triangles 5414 a[0]= (Real8) dete[0]/ tb.det; 5415 a[1]= (Real8) dete[1] / tb.det; 5416 a[2] = (Real8) dete[2] / tb.det; 5417 inside = 1; 5418 return Th.Number(tb); 5419 } 5420 else 5421 { 5422 inside = 0; 5423 double aa,bb; 5424 TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb); 5425 int k = ta; 5426 Triangle * tc = ta; 5427 if (!tc->link) 5428 { ta = ta.Adj(); 5429 tc=ta; 5430 k = ta; 5431 Exchange(aa,bb); 5432 if (!tc->link){ 5433 throw ErrorException(__FUNCT__,exprintf("!tc->link")); 5434 } 5435 } 5436 a[VerticesOfTriangularEdge[k][0]] = aa; 5437 a[VerticesOfTriangularEdge[k][1]] = bb; 5438 a[OppositeVertex[k]] = 1- aa -bb; 5439 return Th.Number(tc); 5440 } 5441 } 5442 /*}}}1*/ 5443 /*FUNCTION CloseBoundaryEdge{{{1*/ 5444 TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) { 5445 int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; 5446 int dir=0; 5447 if (k<0){ 5448 throw ErrorException(__FUNCT__,exprintf("k<0")); 5449 } 5450 int kkk=0; 5451 Icoor2 IJ_IA,IJ_AJ; 5452 TriangleAdjacent edge(t,OppositeEdge[k]); 5453 for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) { 5454 kkk++; 5455 if (kkk>=1000){ 5456 throw ErrorException(__FUNCT__,exprintf("kkk>=1000")); 5457 } 5458 Vertex &vI = *edge.EdgeVertex(0); 5459 Vertex &vJ = *edge.EdgeVertex(1); 5460 I2 I=vI, J=vJ, IJ= J-I; 5461 IJ_IA = (IJ ,(A-I)); 5462 if (IJ_IA<0) { 5463 if (dir>0) {a=1;b=0;return edge;}// change of signe => I 5464 else {dir=-1; 5465 continue;}};// go in direction i 5466 IJ_AJ = (IJ ,(J-A)); 5467 if (IJ_AJ<0) { 5468 if(dir<0) {a=0;b=1;return edge;} 5469 else {dir = 1; 5470 continue;}}// go in direction j 5471 double IJ2 = IJ_IA + IJ_AJ; 5472 if (IJ2==0){ 5473 throw ErrorException(__FUNCT__,exprintf("IJ2==0")); 5474 } 5475 a= IJ_AJ/IJ2; 5476 b= IJ_IA/IJ2; 5477 return edge; 5478 } 5479 } 5480 /*}}}1*/ 5481 /*FUNCTION CloseBoundaryEdgeV2{{{1*/ 5482 TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b) { 5483 // walk around the vertex 5484 // version 2 for remove the probleme if we fill the hole 5485 //int bug=1; 5486 // Triangle *torigine = t; 5487 // restart: 5488 // int dir=0; 5489 if (t->link != 0){ 5490 throw ErrorException(__FUNCT__,exprintf("t->link != 0")); 5491 } 5492 // to have a starting edges 5493 // try the 3 edge bourna-- in case of internal hole 5494 // and choice the best 5495 // 5496 // the probleme is in case of the fine and long internal hole 5497 // for exemple neart the training edge of a wing 5498 Vertex * s=0,*s1=0, *s0=0; 5499 Icoor2 imax = MaxICoor22; 5500 Icoor2 l0 = imax,l1 = imax; 5501 double dd2 = imax;// infinity 5502 TriangleAdjacent er; 5503 int cas=-2; 5504 for (int j=0;j<3;j++) 5505 { 5506 TriangleAdjacent ta=t->FindBoundaryEdge(j); 5507 if (! (Triangle *) ta) continue; 5508 s0 = ta.EdgeVertex(0); 5509 s1 = ta.EdgeVertex(1); 5510 I2 A = * s0; 5511 I2 B = *ta.EdgeVertex(1); 5512 I2 AB = B-A,AC=C-A,BC=B-C; 5513 Icoor2 ACAC = (AC,AC), BCBC = (BC,BC); 5514 Icoor2 AB2 = Norme2_2(AB); // ||AB||^2 5515 Icoor2 ABAC = (AB,AC); // AB.AC| 5516 5517 double d2; 5518 if ( ABAC < 0 ) // DIST A 5519 { 5520 if ( (d2=(double) ACAC) < dd2) 5521 { 5522 er = ta; 5523 l0 = ACAC; 5524 l1 = BCBC; 5525 cas = 0; 5526 s = s0; 5527 } 5528 } 5529 else if (ABAC > AB2) // DIST B 5530 { 5531 if ( (d2=(double) BCBC) < dd2) 5532 { 5533 dd2 = d2; 5534 er = Adj(ta); // other direction 5535 l0 = BCBC; 5536 l1 = ACAC; 5537 cas = 1; 5538 s = s1; 5539 } 5540 } 5541 else // DIST AB 5542 { 5543 5544 double det_2 = (double) Det(AB,AC); 5545 det_2 *= det_2; // square of area*2 of triangle ABC 5546 d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC 5547 5548 if (d2 < dd2) 5549 { 5550 dd2 = d2; 5551 er = ta; 5552 l0 = (AC,AC); 5553 l1 = (BC,BC); 5554 s = 0; 5555 cas = -1; 5556 b = ((double) ABAC/(double) AB2); 5557 a = 1 - b; 5558 } 5559 } 5560 } 5561 if (cas ==-2){ 5562 throw ErrorException(__FUNCT__,exprintf("cas==-2")); 5563 } 5564 // l1 = ||C s1|| , l0 = ||C s0|| 5565 // where s0,s1 are the vertex of the edge er 5566 5567 if ( s) 5568 { 5569 t=er; 5570 TriangleAdjacent edge(er); 5571 5572 int kkk=0; 5573 int linkp = t->link == 0; 5574 5575 Triangle * tt=t=edge=Adj(Previous(edge)); 5576 do { // loop over vertex s 5577 kkk++; 5578 if (edge.EdgeVertex(0)!=s && kkk>=10000){ 5579 throw ErrorException(__FUNCT__,exprintf("edge.EdgeVertex(0)!=s && kkk>=10000")); 5580 } 5581 5582 int link = tt->link == 0; 5583 if ((link + linkp) == 1) 5584 { // a boundary edge 5585 Vertex * st = edge.EdgeVertex(1); 5586 I2 I=*st; 5587 Icoor2 ll = Norme2_2 (C-I); 5588 if (ll < l1) { // the other vertex is neart 5589 s1=st; 5590 l1=ll; 5591 er = edge; 5592 if(ll<l0) { // change of direction -- 5593 s1=s; 5594 l1=l0; 5595 s=st; 5596 l0=ll; 5597 t=tt; 5598 edge=Adj(edge); 5599 link=linkp; 5600 er = edge; 5601 } 5602 } 5603 } 5604 5605 linkp=link; 5606 edge=Adj(Previous(edge)); 5607 tt = edge; 5608 } while (t!=tt); 5609 5610 if (!(Triangle *) er){ 5611 throw ErrorException(__FUNCT__,exprintf("!(Triangle *) er")); 5612 } 5613 I2 A((I2)*er.EdgeVertex(0)); 5614 I2 B((I2)*er.EdgeVertex(1)); 5615 I2 AB=B-A,AC=C-A,CB=B-C; 5616 double aa = (double) (AB,AC); 5617 double bb = (double) (AB,CB); 5618 if (aa<0) a=1,b=0; 5619 else if(bb<0) a=0,b=1; 5620 else 5621 { 5622 a = bb/(aa+bb); 5623 b = aa/(aa+bb); 5624 } 5625 } 5626 return er; 5627 } 5628 /*}}}1*/ 5629 /*FUNCTION AGoodNumberPrimeWith{{{1*/ 5630 Int4 AGoodNumberPrimeWith(Int4 n){ 5631 5632 //list of big prime numbers 5633 const Int4 BigPrimeNumber[] ={ 567890359L, 5634 567890431L, 567890437L, 567890461L, 567890471L, 5635 567890483L, 567890489L, 567890497L, 567890507L, 5636 567890591L, 567890599L, 567890621L, 567890629L , 0}; 5637 5638 //initialize o and pi 5639 Int4 o =0; 5640 Int4 pi=BigPrimeNumber[1]; 5641 5642 //loop until BigPrimeNumber[i]==0 (end of BigPrimeNumber) 5643 for (int i=0; BigPrimeNumber[i]; i++){ 5644 5645 //compute r, rest of the remainder of the division of BigPrimeNumber[i] by n 5646 Int4 r = BigPrimeNumber[i] % n; 5647 5648 /*compute oo = min ( r , n-r , |n - 2r|, |n-3r|)*/ 5649 Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r))); 5650 if ( o < oo){ 5651 o=oo; 5652 pi=BigPrimeNumber[i]; 5653 } 5654 } 5655 return pi; 5656 } 5657 /*}}}1*/ 5658 /*FUNCTION SwapForForcingEdge{{{1*/ 5659 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) { 5660 // l'arete ta coupe l'arete pva pvb 5661 // de cas apres le swap sa coupe toujours 5662 // on cherche l'arete suivante 5663 // on suppose que detsa >0 et detsb <0 5664 // attention la routine echange pva et pvb 5665 5666 if(tt1.Locked()) return 0; // frontiere croise 5667 5668 TriangleAdjacent tt2 = Adj(tt1); 5669 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 5670 Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 5671 if ( a1<0 || a1>=3 ){ 5672 throw ErrorException(__FUNCT__,exprintf("a1<0 || a1>=3")); 5673 } 5674 5675 Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 5676 Vertex & s1= (*t1)[OppositeVertex[a1]]; 5677 Vertex & s2= (*t2)[OppositeVertex[a2]]; 5678 5679 5680 Icoor2 dets2 = det(*pva,*pvb,s2); 5681 Icoor2 det1=t1->det , det2=t2->det ; 5682 Icoor2 detT = det1+det2; 5683 if ((det1<=0 ) || (det2<=0)){ 5684 throw ErrorException(__FUNCT__,exprintf("(det1<=0 ) || (det2<=0)")); 5685 } 5686 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb 5687 throw ErrorException(__FUNCT__,exprintf("(detsa>=0) || (detsb<=0)")); 5688 } 5689 Icoor2 ndet1 = bamg::det(s1,sa,s2); 5690 Icoor2 ndet2 = detT - ndet1; 5691 5692 int ToSwap =0; //pas de swap 5693 if ((ndet1 >0) && (ndet2 >0)) 5694 { // on peut swaper 5695 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 5696 ToSwap =1; 5697 else // swap alleatoire 5698 if (BinaryRand()) 5699 ToSwap =2; 5700 } 5701 if (ToSwap) NbSwap++, 5702 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 5703 5704 int ret=1; 5705 5706 if (dets2 < 0) {// haut 5707 dets1 = ToSwap ? dets1 : detsa ; 5708 detsa = dets2; 5709 tt1 = Previous(tt2) ;} 5710 else if (dets2 > 0){// bas 5711 dets1 = ToSwap ? dets1 : detsb ; 5712 detsb = dets2; 5713 //xxxx tt1 = ToSwap ? tt1 : Next(tt2); 5714 if(!ToSwap) tt1 = Next(tt2); 5715 } 5716 else { // changement de sens 5717 ret = -1; 5718 Exchange(pva,pvb); 5719 Exchange(detsa,detsb); 5720 Exchange(dets1,dets2); 5721 Exchange(tt1,tt2); 5722 dets1=-dets1; 5723 dets2=-dets2; 5724 detsa=-detsa; 5725 detsb=-detsb; 5726 5727 if (ToSwap) 5728 if (dets2 < 0) {// haut 5729 dets1 = (ToSwap ? dets1 : detsa) ; 5730 detsa = dets2; 5731 tt1 = Previous(tt2) ;} 5732 else if (dets2 > 0){// bas 5733 dets1 = (ToSwap ? dets1 : detsb) ; 5734 detsb = dets2; 5735 if(!ToSwap) tt1 = Next(tt2); 5736 } 5737 else {// on a fin ??? 5738 tt1 = Next(tt2); 5739 ret =0;} 5740 5741 } 5742 return ret; 5743 } 5744 /*}}}1*/ 5745 /*FUNCTION ForceEdge{{{1*/ 5746 5747 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) { 5748 int NbSwap =0; 5749 if (!a.t || !b.t){ // the 2 vertex is in a mesh 5750 throw ErrorException(__FUNCT__,exprintf("!a.t || !b.t")); 5751 } 5752 int k=0; 5753 taret=TriangleAdjacent(0,0); // erreur 5754 5755 TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); 5756 Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 5757 // we turn around a in the direct sens 5758 5759 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; 5760 if(v2) // normal case 5761 det2 = det(*v2,a,b); 5762 else { // no chance infini vertex try the next 5763 tta= Previous(Adj(tta)); 5764 v2 = tta.EdgeVertex(0); 5765 vbegin =v2; 5766 if (!v2){ 5767 throw ErrorException(__FUNCT__,exprintf("!v2")); 5768 } 5769 det2 = det(*v2,a,b); 5770 } 5771 5772 while (v2 != &b) { 5773 TriangleAdjacent tc = Previous(Adj(tta)); 5774 v1 = v2; 5775 v2 = tc.EdgeVertex(0); 5776 det1 = det2; 5777 det2 = v2 ? det(*v2,a,b): det2; 5778 5779 if((det1 < 0) && (det2 >0)) { 5780 // try to force the edge 5781 Vertex * va = &a, *vb = &b; 5782 tc = Previous(tc); 5783 if (!v1 || !v2){ 5784 throw ErrorException(__FUNCT__,exprintf("!v1 || !v2")); 5785 } 5786 Icoor2 detss = 0,l=0,ks; 5787 while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 5788 if(l++ > 10000000) { 5789 throw ErrorException(__FUNCT__,exprintf("Loop in forcing Egde, nb de swap=%i, nb of try swap (%i) too big",NbSwap,l)); 5790 } 5791 Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 5792 if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { 5793 tc.SetLock(); 5794 a.Optim(1,0); 5795 b.Optim(1,0); 5796 taret = tc; 5797 return NbSwap; 5798 } 5799 else 5800 { 5801 taret = tc; 5802 return -2; // error boundary is crossing 5803 } 5804 } 5805 tta = tc; 5806 k++; 5807 if (k>=2000){ 5808 throw ErrorException(__FUNCT__,exprintf("k>=2000")); 5809 } 5810 if ( vbegin == v2 ) return -1;// error 5811 } 5812 5813 tta.SetLock(); 5814 taret=tta; 5815 a.Optim(1,0); 5816 b.Optim(1,0); 5817 return NbSwap; 5818 } 5819 /*}}}1*/ 5820 5361 5821 } -
issm/trunk/src/c/Bamgx/shared/shared.h
r3237 r3239 7 7 #define _SHAREDBamg_H_ 8 8 9 #include "AGoodNumberPrimeWith.h"10 #include "CloseBoundaryEdge.h"11 #include "CloseBoundaryEdgeV2.h"12 #include "ForceEdge.h"13 #include "shared.h"14 #include "SwapForForcingEdge.h"15 #include "swap.h"16 9 #include "TheVertex.h" 17 10 #include "FindTriangleAdjacent.h"
Note:
See TracChangeset
for help on using the changeset viewer.
