![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3238
- Timestamp:
- 03/09/10 16:44:03 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 3 edited
-
BamgObjects.h (modified) (1 diff)
-
meshtype.h (modified) (2 diffs)
-
objects/Triangles.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/BamgObjects.h
r3237 r3238 34 34 35 35 /*INLINE functions{{{1*/ 36 // to sort in descending order37 template<class T> inline void HeapSort(T *c,long n){38 /*Intermediary*/39 int l,j,r,i;40 T crit;41 c--; //the array must starts at 1 and not 042 if(n<=1) return; //return if size <=143 l=n/2+1; //initialize l and r44 r=n;45 for(;;){46 if(l<=1){47 crit =c[r];48 c[r--]=c[1];49 if (r==1){c[1]=crit; return;}50 }51 else crit = c[--l];52 j=l;53 for(;;){54 i=j;55 j=2*j;56 if (j>r) {c[i]=crit;break;}57 if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value)58 if (crit < c[j]) c[i]=c[j]; //c[j] > crit -> stock this large value in i(<j)59 else{c[i]=crit;break;} //c[j] < crit -> stock crit in i (<j)60 }61 }62 }63 // to sort in descending order and return ordering64 template<class T> inline void HeapSort(int** porder,T* c,int n){65 /*Intermediary*/66 int l,j,r,i;67 T crit;68 int pos;69 int* order=NULL;70 order=(int*)xmalloc(n*sizeof(int));71 for(i=0;i<n;i++) order[i]=i+1;72 c--; //the array must starts at 1 and not 073 order--;74 if(n<=1) return; //return if size <=175 l=n/2+1; //initialize l and r76 r=n;77 for(;;){78 if(l<=1){79 crit =c[r]; pos=order[r];80 c[r--]=c[1]; order[r+1]=order[1];81 if (r==1){82 c[1]=crit; order[1]=pos;83 order++;84 *porder=order;85 return;86 }87 }88 else {crit=c[--l]; pos=order[l];}89 j=l;90 for(;;){91 i=j;92 j=2*j;93 if (j>r) {c[i]=crit;order[i]=pos;break;}94 if ((j<r) && (c[j] < c[j+1]))j++;95 if (crit < c[j]) {c[i]=c[j];order[i]=order[j];}96 else{c[i]=crit;order[i]=pos;break;}97 }98 }99 }100 36 inline Real8 det3x3(Real8 A[3] ,Real8 B[3],Real8 C[3]){ 101 37 return A[0]*( B[1]*C[2]-B[2]*C[1]) -
issm/trunk/src/c/Bamgx/meshtype.h
r3234 r3238 58 58 template<class T> inline T Max3 (const T &a,const T & b,const T & c){return Max(Max(a,b),c);} 59 59 template<class T> inline T Min3 (const T &a,const T & b,const T & c){return Min(Min(a,b),c);} 60 template<class T> inline void HeapSort(T *c,long n){ 61 /*Intermediary*/ 62 int l,j,r,i; 63 T crit; 64 c--; //the array must starts at 1 and not 0 65 if(n<=1) return; //return if size <=1 66 l=n/2+1; //initialize l and r 67 r=n; 68 for(;;){ 69 if(l<=1){ 70 crit =c[r]; 71 c[r--]=c[1]; 72 if (r==1){c[1]=crit; return;} 73 } 74 else crit = c[--l]; 75 j=l; 76 for(;;){ 77 i=j; 78 j=2*j; 79 if (j>r) {c[i]=crit;break;} 80 if ((j<r) && (c[j] < c[j+1])) j++;//c[j+1]> c[j] -> take j+1 instead of j (larger value) 81 if (crit < c[j]) c[i]=c[j]; //c[j] > crit -> stock this large value in i(<j) 82 else{c[i]=crit;break;} //c[j] < crit -> stock crit in i (<j) 83 } 84 } 85 } 86 template<class T> inline void HeapSort(int** porder,T* c,int n){ 87 /*Intermediary*/ 88 int l,j,r,i; 89 T crit; 90 int pos; 91 int* order=NULL; 92 order=(int*)xmalloc(n*sizeof(int)); 93 for(i=0;i<n;i++) order[i]=i+1; 94 c--; //the array must starts at 1 and not 0 95 order--; 96 if(n<=1) return; //return if size <=1 97 l=n/2+1; //initialize l and r 98 r=n; 99 for(;;){ 100 if(l<=1){ 101 crit =c[r]; pos=order[r]; 102 c[r--]=c[1]; order[r+1]=order[1]; 103 if (r==1){ 104 c[1]=crit; order[1]=pos; 105 order++; 106 *porder=order; 107 return; 108 } 109 } 110 else {crit=c[--l]; pos=order[l];} 111 j=l; 112 for(;;){ 113 i=j; 114 j=2*j; 115 if (j>r) {c[i]=crit;order[i]=pos;break;} 116 if ((j<r) && (c[j] < c[j+1]))j++; 117 if (crit < c[j]) {c[i]=c[j];order[i]=order[j];} 118 else{c[i]=crit;order[i]=pos;break;} 119 } 120 } 121 } 60 122 61 123 //Inline functions … … 77 139 } 78 140 79 80 141 } 81 142 #endif -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r3234 r3238 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 triangles5366 //5367 // sb sb5368 // / | \ / \ !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 sa5377 // -------------------------------------------------------------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 sb5383 (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa5384 // mise a jour des 2 adjacences externes5385 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 gauche5390 taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap());5391 // externe bas droite5392 taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap());5393 // remove the Mark UnMarkSwap5394 t1->SetUnMarkUnSwap(ap1);5395 t2->SetUnMarkUnSwap(ap2);5396 // interne5397 tas1.SetAdj2(tas2);5398 5399 t1->det = det1;5400 t2->det = det2;5401 5402 t1->SetTriangleContainingTheVertex();5403 t2->SetTriangleContainingTheVertex();5404 } // end swap5405 /*}}}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 triangles5414 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 else5421 {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 => I5464 else {dir=-1;5465 continue;}};// go in direction i5466 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 j5471 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 vertex5484 // version 2 for remove the probleme if we fill the hole5485 //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 edges5493 // try the 3 edge bourna-- in case of internal hole5494 // and choice the best5495 //5496 // the probleme is in case of the fine and long internal hole5497 // for exemple neart the training edge of a wing5498 Vertex * s=0,*s1=0, *s0=0;5499 Icoor2 imax = MaxICoor22;5500 Icoor2 l0 = imax,l1 = imax;5501 double dd2 = imax;// infinity5502 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||^25515 Icoor2 ABAC = (AB,AC); // AB.AC|5516 5517 double d2;5518 if ( ABAC < 0 ) // DIST A5519 {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 B5530 {5531 if ( (d2=(double) BCBC) < dd2)5532 {5533 dd2 = d2;5534 er = Adj(ta); // other direction5535 l0 = BCBC;5536 l1 = ACAC;5537 cas = 1;5538 s = s1;5539 }5540 }5541 else // DIST AB5542 {5543 5544 double det_2 = (double) Det(AB,AC);5545 det_2 *= det_2; // square of area*2 of triangle ABC5546 d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC5547 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 er5566 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 s5577 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 edge5585 Vertex * st = edge.EdgeVertex(1);5586 I2 I=*st;5587 Icoor2 ll = Norme2_2 (C-I);5588 if (ll < l1) { // the other vertex is neart5589 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 else5621 {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 numbers5633 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 pi5639 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 n5646 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 pvb5661 // de cas apres le swap sa coupe toujours5662 // on cherche l'arete suivante5663 // on suppose que detsa >0 et detsb <05664 // attention la routine echange pva et pvb5665 5666 if(tt1.Locked()) return 0; // frontiere croise5667 5668 TriangleAdjacent tt2 = Adj(tt1);5669 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent5670 Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles5671 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,bb5687 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 swap5693 if ((ndet1 >0) && (ndet2 >0))5694 { // on peut swaper5695 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))5696 ToSwap =1;5697 else // swap alleatoire5698 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) {// haut5707 dets1 = ToSwap ? dets1 : detsa ;5708 detsa = dets2;5709 tt1 = Previous(tt2) ;}5710 else if (dets2 > 0){// bas5711 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 sens5717 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) {// haut5729 dets1 = (ToSwap ? dets1 : detsa) ;5730 detsa = dets2;5731 tt1 = Previous(tt2) ;}5732 else if (dets2 > 0){// bas5733 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 mesh5750 throw ErrorException(__FUNCT__,exprintf("!a.t || !b.t"));5751 }5752 int k=0;5753 taret=TriangleAdjacent(0,0); // erreur5754 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 sens5758 5759 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1;5760 if(v2) // normal case5761 det2 = det(*v2,a,b);5762 else { // no chance infini vertex try the next5763 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 edge5781 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 else5800 {5801 taret = tc;5802 return -2; // error boundary is crossing5803 }5804 }5805 tta = tc;5806 k++;5807 if (k>=2000){5808 throw ErrorException(__FUNCT__,exprintf("k>=2000"));5809 }5810 if ( vbegin == v2 ) return -1;// error5811 }5812 5813 tta.SetLock();5814 taret=tta;5815 a.Optim(1,0);5816 b.Optim(1,0);5817 return NbSwap;5818 }5819 /*}}}1*/5820 5821 5361 }
Note:
See TracChangeset
for help on using the changeset viewer.
