![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3326
- Timestamp:
- 03/24/10 09:54:58 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx/objects
- Files:
-
- 4 edited
-
GeometricalEdge.h (modified) (4 diffs)
-
Geometry.cpp (modified) (1 diff)
-
Triangles.cpp (modified) (8 diffs)
-
Triangles.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/objects/GeometricalEdge.h
r3280 r3326 28 28 int DirAdj[2]; 29 29 int flag ; 30 GeometricalEdge* link; // if Cracked() or Equi()31 30 32 31 //Operators … … 40 39 int Tg(int i) const{return i==0 ? TgA() : TgB();} 41 40 int Cracked() const {return flag & 1; } 42 int Dup() const {return flag & 32; }43 int Equi()const {return flag & 2; }44 int ReverseEqui()const {return flag & 128;}45 41 int TgA()const {return flag & 4; } 46 42 int TgB()const {return flag & 8; } … … 48 44 int Required() {return flag & 64; } 49 45 void SetCracked() { flag |= 1; } 50 void SetDup() { flag |= 32; } // not a real edge51 void SetEqui() { flag |= 2; }52 46 void SetTgA() { flag |=4; } 53 47 void SetTgB() { flag |=8; } … … 55 49 void SetUnMark() { flag &= 1007 /* 1023-16*/;} 56 50 void SetRequired() { flag |= 64; } 57 void SetReverseEqui() { flag |= 128;}58 51 void Set(const GeometricalEdge & rec,const Geometry & Th ,Geometry & ThNew); 59 52 }; -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r3325 r3326 739 739 } 740 740 for(int i=0;i<NbOfCurves ;i++){ 741 GeometricalEdge * be=curves[i].be, *eqbe=be ->link;742 //GeometricalEdge * ee=curves[i].ee, *eqee=be ->link;741 GeometricalEdge * be=curves[i].be, *eqbe=be; 742 //GeometricalEdge * ee=curves[i].ee, *eqee=be; 743 743 curves[i].master=true; 744 if(be->Equi() || be->ReverseEqui() ){745 if (!eqbe){746 throw ErrorException(__FUNCT__,exprintf("!eqbe"));747 }748 int nc = eqbe->CurveNumber;749 if (i==nc){750 throw ErrorException(__FUNCT__,exprintf("i==nc"));751 }752 curves[i].next=curves[nc].next;753 curves[i].master=false;754 curves[nc].next=curves+i;755 if(be->ReverseEqui())756 curves[i].Reverse();757 }758 744 } 759 745 -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r3325 r3326 379 379 } 380 380 } 381 else{382 if(verbose>5) printf(" no Quadrilaterals found\n");383 }384 381 385 382 //VerticesOnGeometricEdge … … 397 394 } 398 395 } 399 else{400 if(verbose>5) printf(" no VertexOnGeometricEdge found\n");401 }402 396 403 397 //VerticesOnGeometricVertex … … 412 406 VerticesOnGeomVertex[i]=VertexOnGeom(vertices[i1],Gh.vertices[i2]); 413 407 } 414 }415 else{416 if(verbose>5) printf(" no VertexOnGeometricVertex found\n");417 408 } 418 409 … … 479 470 } 480 471 } 481 else{482 if(verbose>5) printf(" no Edges found\n");483 }484 472 485 473 //EdgeOnGeometricEdge … … 498 486 } 499 487 } 500 else{501 if(verbose>5) printf(" no EdgesOnGeometricEdge found\n");502 }503 488 504 489 //hVertices … … 510 495 } 511 496 } 512 }513 else{514 if(verbose>5) printf(" no hVertices found\n");515 497 } 516 498 … … 529 511 subdomains[i].head = triangles+head; 530 512 } 531 }532 else{533 if(verbose>5) printf(" no SubDomains found\n");534 513 } 535 514 … … 2655 2634 GeometricalEdge &ei=Gh.edges[i]; 2656 2635 2657 // a good curve (not a duplicate) 2658 if (!ei.Dup()){ 2659 2660 for(int j=0;j<2;j++) { 2661 2662 /*The first time, no edge is marked but this might change during the loop*/ 2663 if (!ei.Mark() && ei[j].Required()){ 2664 2665 long nbvend=0; 2666 Edge* PreviousNewEdge=NULL; 2667 2668 lstep = -1; 2669 2670 /*If Edge is required*/ 2671 if(ei.Required()){ 2672 //do not create internal points if required (take it as is) 2673 if (j==0){ 2674 if(step==0) nbe++; 2675 else{ 2676 e=&ei; 2677 a=ei(0)->The(); 2678 b=ei(1)->The(); 2679 2680 //check that edges has been allocated 2681 if (!edges){ 2682 throw ErrorException(__FUNCT__,exprintf("edges has not been allocated...")); 2683 } 2684 edges[nbe].v[0]=a->to; 2685 edges[nbe].v[1]=b->to;; 2686 edges[nbe].ref = e->ref; 2687 edges[nbe].onGeometry = e; 2688 edges[nbe].adj[0] = 0; 2689 edges[nbe].adj[1] = 0; 2690 nbe++; 2636 for(int j=0;j<2;j++) { 2637 2638 /*The first time, no edge is marked but this might change during the loop*/ 2639 if (!ei.Mark() && ei[j].Required()){ 2640 2641 long nbvend=0; 2642 Edge* PreviousNewEdge=NULL; 2643 2644 lstep = -1; 2645 2646 /*If Edge is required*/ 2647 if(ei.Required()){ 2648 //do not create internal points if required (take it as is) 2649 if (j==0){ 2650 if(step==0) nbe++; 2651 else{ 2652 e=&ei; 2653 a=ei(0)->The(); 2654 b=ei(1)->The(); 2655 2656 //check that edges has been allocated 2657 if (!edges){ 2658 throw ErrorException(__FUNCT__,exprintf("edges has not been allocated...")); 2691 2659 } 2660 edges[nbe].v[0]=a->to; 2661 edges[nbe].v[1]=b->to;; 2662 edges[nbe].ref = e->ref; 2663 edges[nbe].onGeometry = e; 2664 edges[nbe].adj[0] = 0; 2665 edges[nbe].adj[1] = 0; 2666 nbe++; 2692 2667 } 2693 2668 } 2694 2695 /*If Edge is not required: on a curve*/ 2696 else { 2697 for ( int kstep=0;kstep<=step;kstep++){ 2698 //step=0, do nothing 2699 //step=1, compute the length of the curve 2700 //step=2 create the points and edge 2701 PreviousNewEdge=0; 2702 NbNewPoints=0; 2703 NbEdgeCurve=0; 2704 if (nbvend>=nbvx){ 2705 throw ErrorException(__FUNCT__,exprintf("nbvend>=nbvx")); 2669 } 2670 2671 /*If Edge is not required: on a curve*/ 2672 else { 2673 for ( int kstep=0;kstep<=step;kstep++){ 2674 //step=0, do nothing 2675 //step=1, compute the length of the curve 2676 //step=2 create the points and edge 2677 PreviousNewEdge=0; 2678 NbNewPoints=0; 2679 NbEdgeCurve=0; 2680 if (nbvend>=nbvx){ 2681 throw ErrorException(__FUNCT__,exprintf("nbvend>=nbvx")); 2682 } 2683 lcurve =0; 2684 s = lstep; 2685 2686 // i = edge number, j=[0;1] vertex number in edge 2687 2688 k=j; // k = vertex number in edge (0 or 1) 2689 e=&ei; // e = reference of current edge 2690 a=ei(k)->The(); // a = pointer toward the kth vertex of the current edge 2691 va = a->to; // va = pointer toward vertex associated 2692 e->SetMark(); // Mark edge 2693 2694 //if SameGeo We have go to the background geometry 2695 //to find the discretisation of the curve 2696 for(;;){ 2697 k = 1-k; 2698 b = (*e)(k)->The();// b = pointer toward the other vertex of the current edge 2699 AB= b->r - a->r; // AB = vector of the current edge 2700 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A 2701 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B 2702 double ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric 2703 2704 /* We are now creating the edges of the mesh from the 2705 * geometrical edge selected above. 2706 * The edge will be divided according to the metric 2707 * previously computed and cannot be divided more 2708 * than 10 times (MaxSubEdge). */ 2709 2710 //By default, there is only one subedge that is the geometrical edge itself 2711 int NbSubEdge = 1; 2712 2713 //initialize lSubEdge, holding the length of each subedge (cannot be higher than 10) 2714 double lSubEdge[MaxSubEdge]; 2715 2716 //Build Subedges according to the edge length 2717 //if ledge < 1.5 (between one and 2), take the edge as is 2718 if (ledge < 1.5) lSubEdge[0] = ledge; 2719 //else, divide the edge 2720 else { 2721 //compute number of subedges (division of the edge) 2722 NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5)); 2723 //A and B are the position of points on the edge 2724 R2 A,B; 2725 A=a->r; 2726 Metric MAs=MA,MBs; 2727 ledge=0; 2728 double x =0, xstep= 1./NbSubEdge; 2729 for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { 2730 x += xstep; 2731 B = e->F(k ? x : 1-x); 2732 MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB); 2733 AB = A-B; 2734 lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2); 2735 } 2706 2736 } 2707 lcurve =0; 2708 s = lstep; 2709 2710 // i = edge number, j=[0;1] vertex number in edge 2711 2712 k=j; // k = vertex number in edge (0 or 1) 2713 e=&ei; // e = reference of current edge 2714 a=ei(k)->The(); // a = pointer toward the kth vertex of the current edge 2715 va = a->to; // va = pointer toward vertex associated 2716 e->SetMark(); // Mark edge 2717 2718 //if SameGeo We have go to the background geometry 2719 //to find the discretisation of the curve 2720 for(;;){ 2721 k = 1-k; 2722 b = (*e)(k)->The();// b = pointer toward the other vertex of the current edge 2723 AB= b->r - a->r; // AB = vector of the current edge 2724 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A 2725 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B 2726 double ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric 2727 2728 /* We are now creating the edges of the mesh from the 2729 * geometrical edge selected above. 2730 * The edge will be divided according to the metric 2731 * previously computed and cannot be divided more 2732 * than 10 times (MaxSubEdge). */ 2733 2734 //By default, there is only one subedge that is the geometrical edge itself 2735 int NbSubEdge = 1; 2736 2737 //initialize lSubEdge, holding the length of each subedge (cannot be higher than 10) 2738 double lSubEdge[MaxSubEdge]; 2739 2740 //Build Subedges according to the edge length 2741 //if ledge < 1.5 (between one and 2), take the edge as is 2742 if (ledge < 1.5) lSubEdge[0] = ledge; 2743 //else, divide the edge 2744 else { 2745 //compute number of subedges (division of the edge) 2746 NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5)); 2747 //A and B are the position of points on the edge 2748 R2 A,B; 2749 A=a->r; 2750 Metric MAs=MA,MBs; 2751 ledge=0; 2752 double x =0, xstep= 1./NbSubEdge; 2753 for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { 2754 x += xstep; 2755 B = e->F(k ? x : 1-x); 2756 MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB); 2757 AB = A-B; 2758 lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2); 2759 } 2737 2738 double lcurveb = lcurve+ ledge ; 2739 2740 //Now, create corresponding points 2741 while (lcurve<=s && s <= lcurveb && nbv < nbvend){ 2742 2743 double ss = s-lcurve; 2744 2745 /*Find the SubEdge containing ss using Dichotomy*/ 2746 2747 int kk0=-1,kk1=NbSubEdge-1,kkk; 2748 double ll0=0,ll1=ledge,llk; 2749 while (kk1-kk0>1){ 2750 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) 2751 kk1=kkk,ll1=llk; 2752 else 2753 kk0=kkk,ll0=llk; 2760 2754 } 2761 2762 double lcurveb = lcurve+ ledge ; 2763 2764 //Now, create corresponding points 2765 while (lcurve<=s && s <= lcurveb && nbv < nbvend){ 2766 2767 double ss = s-lcurve; 2768 2769 /*Find the SubEdge containing ss using Dichotomy*/ 2770 2771 int kk0=-1,kk1=NbSubEdge-1,kkk; 2772 double ll0=0,ll1=ledge,llk; 2773 while (kk1-kk0>1){ 2774 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) 2775 kk1=kkk,ll1=llk; 2776 else 2777 kk0=kkk,ll0=llk; 2778 } 2779 if (kk1 == kk0){ 2780 throw ErrorException(__FUNCT__,exprintf("kk1 == kk0")); 2781 } 2782 2783 double sbb = (ss-ll0 )/(ll1-ll0); 2784 double bb = (kk1+sbb)/NbSubEdge, aa=1-bb; 2785 2786 // new vertex on edge 2787 vb = &vertices[nbv++]; 2788 vb->m = Metric(aa,a->m,bb,b->m); 2789 vb->ReferenceNumber = e->ref; 2790 vb->DirOfSearch =NoDirOfSearch; 2791 double abcisse = k ? bb : aa; 2792 vb->r = e->F( abcisse ); 2793 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); 2794 2795 // to take into account the direction of the edge 2796 s += lstep; 2797 edges[nbe].v[0]=va; 2798 edges[nbe].v[1]=vb; 2799 edges[nbe].ref =e->ref; 2800 edges[nbe].onGeometry = e; 2801 edges[nbe].adj[0] = PreviousNewEdge; 2802 if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe]; 2803 PreviousNewEdge=edges+nbe; 2804 nbe++; 2805 va = vb; 2755 if (kk1 == kk0){ 2756 throw ErrorException(__FUNCT__,exprintf("kk1 == kk0")); 2806 2757 } 2807 lcurve = lcurveb; 2808 e->SetMark(); 2809 a=b; 2810 if (b->Required() ) break; 2811 int kprev=k; 2812 k = e->DirAdj[kprev];// next vertices 2813 e = e->Adj[kprev]; 2814 if (!e){ 2815 throw ErrorException(__FUNCT__,exprintf("!e")); 2816 } 2817 }// for(;;) 2818 vb = b->to; 2819 NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1); 2820 NbNewPoints = NbEdgeCurve-1; 2821 if(!kstep){ 2822 NbVerticesOnGeomEdge0 += NbNewPoints; 2823 NbOfCurves++; 2758 2759 double sbb = (ss-ll0 )/(ll1-ll0); 2760 double bb = (kk1+sbb)/NbSubEdge, aa=1-bb; 2761 2762 // new vertex on edge 2763 vb = &vertices[nbv++]; 2764 vb->m = Metric(aa,a->m,bb,b->m); 2765 vb->ReferenceNumber = e->ref; 2766 vb->DirOfSearch =NoDirOfSearch; 2767 double abcisse = k ? bb : aa; 2768 vb->r = e->F( abcisse ); 2769 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); 2770 2771 // to take into account the direction of the edge 2772 s += lstep; 2773 edges[nbe].v[0]=va; 2774 edges[nbe].v[1]=vb; 2775 edges[nbe].ref =e->ref; 2776 edges[nbe].onGeometry = e; 2777 edges[nbe].adj[0] = PreviousNewEdge; 2778 if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe]; 2779 PreviousNewEdge=edges+nbe; 2780 nbe++; 2781 va = vb; 2824 2782 } 2825 nbvend=nbv+NbNewPoints; 2826 lstep = lcurve / NbEdgeCurve; 2827 }// end of curve -- 2828 if (edges) { // last edges of the curves 2829 edges[nbe].v[0]=va; 2830 edges[nbe].v[1]=vb; 2831 edges[nbe].ref = e->ref; 2832 edges[nbe].onGeometry = e; 2833 edges[nbe].adj[0] = PreviousNewEdge; 2834 edges[nbe].adj[1] = 0; 2835 if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe]; 2836 nbe++; 2783 lcurve = lcurveb; 2784 e->SetMark(); 2785 a=b; 2786 if (b->Required() ) break; 2787 int kprev=k; 2788 k = e->DirAdj[kprev];// next vertices 2789 e = e->Adj[kprev]; 2790 if (!e){ 2791 throw ErrorException(__FUNCT__,exprintf("!e")); 2792 } 2793 }// for(;;) 2794 vb = b->to; 2795 NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1); 2796 NbNewPoints = NbEdgeCurve-1; 2797 if(!kstep){ 2798 NbVerticesOnGeomEdge0 += NbNewPoints; 2799 NbOfCurves++; 2837 2800 } 2838 else nbe += NbEdgeCurve; 2839 } // end on curve --- 2840 } 2801 nbvend=nbv+NbNewPoints; 2802 lstep = lcurve / NbEdgeCurve; 2803 }// end of curve -- 2804 if (edges) { // last edges of the curves 2805 edges[nbe].v[0]=va; 2806 edges[nbe].v[1]=vb; 2807 edges[nbe].ref = e->ref; 2808 edges[nbe].onGeometry = e; 2809 edges[nbe].adj[0] = PreviousNewEdge; 2810 edges[nbe].adj[1] = 0; 2811 if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe]; 2812 nbe++; 2813 } 2814 else nbe += NbEdgeCurve; 2815 } // end on curve --- 2841 2816 } 2842 2817 } -
issm/trunk/src/c/Bamgx/objects/Triangles.h
r3325 r3326 148 148 //Inline methods 149 149 inline void ReMakeTriangleContainingTheVertex(){ 150 for (int i=0;i<nbv;i++) vertices[i].vint=0, vertices[i].t= 0;150 for (int i=0;i<nbv;i++) vertices[i].vint=0, vertices[i].t=NULL; 151 151 for (int i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex(); 152 152 } … … 162 162 } 163 163 inline void SetVertexFieldOnBTh(){ 164 for (int i=0;i<nbv;i++) vertices[i].onGeometry= 0;164 for (int i=0;i<nbv;i++) vertices[i].onGeometry=NULL; 165 165 for (int j=0;j<NbVertexOnBThVertex;j++) VertexOnBThVertex[j].SetOnBTh(); 166 166 for (int k=0;k<NbVertexOnBThEdge;k++ ) VertexOnBThEdge[k].SetOnBTh();
Note:
See TracChangeset
for help on using the changeset viewer.
