![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 5397
- Timestamp:
- 08/19/10 09:10:32 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Bamg
- Files:
-
- 8 edited
-
Curve.cpp (modified) (1 diff)
-
Curve.h (modified) (1 diff)
-
GeometricalEdge.cpp (modified) (1 diff)
-
Geometry.cpp (modified) (11 diffs)
-
Geometry.h (modified) (1 diff)
-
MatVVP2x2.cpp (modified) (1 diff)
-
Mesh.cpp (modified) (5 diffs)
-
Metric.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Bamg/Curve.cpp
r5149 r5397 11 11 /*Constructors/Destructors*/ 12 12 /*FUNCTION Curve::Curve(){{{1*/ 13 Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0) ,master(true){13 Curve::Curve():be(0),ee(0),kb(0),ke(0),next(0){ 14 14 15 15 } -
issm/trunk/src/c/objects/Bamg/Curve.h
r5130 r5397 16 16 int kb,ke; // begin vetex and end vertex 17 17 Curve* next; // next curve equi to this 18 bool master; // true => of equi curve point on this curve19 18 20 19 //Methods -
issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp
r5151 r5397 152 152 /*FUNCTION GeometricalEdge::SetCracked{{{1*/ 153 153 void GeometricalEdge::SetCracked() { 154 type |= 1; 154 type |= 1;/*=>1st digit to 1*/ 155 155 }/*}}}*/ 156 156 /*FUNCTION GeometricalEdge::SetTgA{{{1*/ 157 157 void GeometricalEdge::SetTgA() { 158 type |=4; 158 type |=4; /*=>2d digit to 1*/ 159 159 }/*}}}*/ 160 160 /*FUNCTION GeometricalEdge::SetTgB{{{1*/ 161 161 void GeometricalEdge::SetTgB() { 162 type |=8; 162 type |=8; /*=> 3d digit to 1*/ 163 163 }/*}}}*/ 164 164 /*FUNCTION GeometricalEdge::SetMark{{{1*/ 165 165 void GeometricalEdge::SetMark() { 166 type |=16; 166 type |=16;/*=> 4th digiy to 1*/ 167 167 }/*}}}*/ 168 168 /*FUNCTION GeometricalEdge::SetUnMark{{{1*/ 169 169 void GeometricalEdge::SetUnMark() { 170 type &= 1007 /* 1023-16 */;170 type &= 1007 /* 1023-16 = 000111110111 => 4th digit to 0*/; 171 171 }/*}}}*/ 172 172 /*FUNCTION GeometricalEdge::SetRequired{{{1*/ 173 173 void GeometricalEdge::SetRequired() { 174 type |= 64; 174 type |= 64;/*=>=>6th digit to 1*/ 175 175 }/*}}}*/ 176 176 /*FUNCTION GeometricalEdge::Tg{{{1*/ -
issm/trunk/src/c/objects/Bamg/Geometry.cpp
r5340 r5397 22 22 Init(); 23 23 ReadGeometry(bamggeom,bamgopts); 24 AfterRead();24 PostRead(); 25 25 } 26 26 /*}}}*/ … … 165 165 } 166 166 167 // definition 167 // definition the default of the given mesh size 168 168 for (i=0;i<nbv;i++) { 169 169 if (vertices[i].color > 0) … … 400 400 401 401 /*Methods*/ 402 /*FUNCTION Geometry:: AfterRead{{{1*/403 void Geometry:: AfterRead(){402 /*FUNCTION Geometry::PostRead{{{1*/ 403 void Geometry::PostRead(){ 404 404 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/ 405 405 … … 606 606 } 607 607 608 / / generation of all the tangents608 /* generation of all the tangents*/ 609 609 for (i=0;i<nbe;i++) { 610 610 R2 AB =edges[i].v[1]->r -edges[i].v[0]->r;// AB = vertex0 -> vertex1 … … 646 646 } 647 647 648 /* generation of all curves (from corner to corner)*/ 649 /*We proceed in 2 steps: first allocate, second build*/ 648 650 for (int step=0;step<2;step++){ 649 651 652 //unmark all edges 650 653 for (i=0;i<nbe;i++) edges[i].SetUnMark(); 651 654 long nb_marked_edges=0; 655 656 //initialize number of curves 652 657 nbcurves = 0; 653 long nbgem=0; 654 655 for (int level=0;level < 2 && nbgem != nbe;level++){ 658 659 for (int level=0;level<2 && nb_marked_edges!=nbe;level++){ 656 660 for (i=0;i<nbe;i++){ 657 GeometricalEdge & ei = edges[i]; 661 662 GeometricalEdge & ei=edges[i]; 658 663 for(j=0;j<2;j++){ 664 /*If current edge ei is unmarked and (level=1 or vertex i is required (corner)): 665 * we do have the first edge of a new curve*/ 659 666 if (!ei.Mark() && (level || ei[j].Required())) { 660 // warning ei.Mark() can be change in loop for(j=0;j<2;j++)661 667 int k0=j,k1; 662 GeometricalEdge *e = &ei;668 GeometricalEdge *e=&ei; 663 669 GeometricalVertex *a=(*e)(k0); // begin 664 670 if(curves){ … … 667 673 } 668 674 int nee=0; 669 for(;;) 675 for(;;){ 670 676 nee++; 671 677 k1 = 1-k0; // next vertex of the edge 672 678 e->SetMark(); 673 nb gem++;679 nb_marked_edges++; 674 680 e->CurveNumber=nbcurves; 675 if(curves) 681 if(curves){ 676 682 curves[nbcurves].ee=e; 677 683 curves[nbcurves].ke=k1; … … 679 685 680 686 GeometricalVertex *b=(*e)(k1); 681 if (a == b || b->Required() ) break; 682 k0 = e->AdjVertexNumber[k1];// vertex in next edge 683 e = e->Adj[k1]; // next edge 687 688 //break if we have reached the other end of the curve 689 if (a==b || b->Required()){ 690 break; 691 } 692 //else: go to next edge (adjacent) 693 else{ 694 k0 = e->AdjVertexNumber[k1];// vertex in next edge 695 e = e->Adj[k1]; // next edge 696 } 684 697 } 685 698 nbcurves++; … … 689 702 } 690 703 } 691 ISSMASSERT(nbgem && nbe); 692 if(step==0){ 693 curves = new Curve[nbcurves]; 694 } 704 ISSMASSERT(nb_marked_edges && nbe); 705 //allocate if first step 706 if(step==0) curves=new Curve[nbcurves]; 695 707 } 696 708 for(int i=0;i<nbcurves ;i++){ 697 709 GeometricalEdge * be=curves[i].be, *eqbe=be; 698 710 //GeometricalEdge * ee=curves[i].ee, *eqee=be; 699 curves[i].master=true;700 711 } 701 712 … … 707 718 } 708 719 /*}}}1*/ 709 /*FUNCTION Geometry::Containing{{{1*/710 GeometricalEdge* Geometry::Containing(const R2 P, GeometricalEdge * start) const {711 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/712 713 GeometricalEdge* on =start,* pon=0;714 // walk with the cos on geometry715 int counter=0;716 while(pon != on){717 counter++;718 ISSMASSERT(counter<100);719 pon = on;720 R2 A= (*on)[0];721 R2 B= (*on)[1];722 R2 AB = B-A;723 R2 AP = P-A;724 R2 BP = P-B;725 if ( (AB,AP) < 0)726 on = on->Adj[0];727 else if ( (AB,BP) > 0)728 on = on->Adj[1];729 else730 return on;731 }732 return on;733 }734 /*}}}1*/735 720 /*FUNCTION Geometry::Echo {{{1*/ 736 737 721 void Geometry::Echo(void){ 738 722 … … 774 758 /*FUNCTION Geometry::MinimalHmin{{{1*/ 775 759 double Geometry::MinimalHmin() { 760 /* coeffIcoor = (2^30-1)/D 761 * We cannot go beyond hmin = D/2^30 because of the quadtree 762 * hmin is therefore approximately 2/coeffIcoor */ 776 763 return 2.0/coefIcoor; 777 764 }/*}}}*/ … … 800 787 return c - curves; 801 788 }/*}}}*/ 789 /*FUNCTION Geometry::Containing{{{1*/ 790 GeometricalEdge* Geometry::Containing(const R2 P, GeometricalEdge * start) const { 791 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/ 792 793 GeometricalEdge* on =start,* pon=0; 794 // walk with the cos on geometry 795 int counter=0; 796 while(pon != on){ 797 counter++; 798 ISSMASSERT(counter<100); 799 pon = on; 800 R2 A= (*on)[0]; 801 R2 B= (*on)[1]; 802 R2 AB = B-A; 803 R2 AP = P-A; 804 R2 BP = P-B; 805 if ( (AB,AP) < 0) 806 on = on->Adj[0]; 807 else if ( (AB,BP) > 0) 808 on = on->Adj[1]; 809 else 810 return on; 811 } 812 return on; 813 } 814 /*}}}1*/ 802 815 /*FUNCTION Geometry::ProjectOnCurve {{{1*/ 803 816 GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const { -
issm/trunk/src/c/objects/Bamg/Geometry.h
r5180 r5397 53 53 void ReadGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts); 54 54 void Init(void); 55 void AfterRead();55 void PostRead(); 56 56 long GetId(const GeometricalVertex &t) const; 57 57 long GetId(const GeometricalVertex *t) const; -
issm/trunk/src/c/objects/Bamg/MatVVP2x2.cpp
r5151 r5397 11 11 /*FUNCTION MatVVP2x2::MatVVP2x2(const Metric M){{{1*/ 12 12 MatVVP2x2::MatVVP2x2(const Metric M){ 13 /* Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MatVVP2x2)*/13 /*From a metric (a11,a21,a22), get eigen values lambda1 and lambda2 and one eigen vector v*/ 14 14 15 /*Intermediaries*/ 15 16 double a11=M.a11,a21=M.a21,a22=M.a22; 16 const double eps = 1.e-5; 17 double c11 = a11*a11, c22 = a22*a22, c21= a21*a21; 18 double b=-a11-a22,c=-c21+a11*a22; 19 double delta = b*b - 4 * c ; 20 double n2=(c11+c22+c21); 17 double norm1,norm2,normM; 18 double delta,b; 21 19 22 //if norm(M)<10^30 -> M=zeros(2,2) 23 if ( n2 < 1e-30) lambda1=lambda2=0,v.x=1,v.y=0; 20 /*To get the eigen values, we must solve the following equation: 21 * | a11 - lambda a21 | 22 * det | | = 0 23 * | a21 a22-lambda | 24 * 25 * We have to solve the following polynom: 26 * lamda^2 + ( -a11 -a22)*lambda + (a11*a22-a21*a21) = 0*/ 24 27 25 //if matrix is already diagonal and has a 0 eigen value 26 else if (delta < eps*n2){ 27 lambda1=lambda2=-b/2, v.x=1,v.y=0; 28 /*Compute polynom determinant*/ 29 b=-a11-a22; 30 delta=b*b - 4*(a11*a22-a21*a21); 31 32 33 /*Compute norm of M to avoid round off errors*/ 34 normM=a11*a11 + a22*a22 + a21*a21; 35 36 /*1: normM too small: eigen values = 0*/ 37 if(normM<1.e-30){ 38 lambda1=0; 39 lambda2=0; 40 v.x=1; 41 v.y=0; 42 } 43 /*2: delta is small -> double root*/ 44 else if (delta < 1.e-5*normM){ 45 lambda1=-b/2; 46 lambda2=-b/2; 47 v.x=1; 48 v.y=0; 49 } 50 /*3: general case -> two roots*/ 51 else{ 52 delta = sqrt(delta); 53 lambda1 = (-b-delta)/2.0; 54 lambda2 = (-b+delta)/2.0; 55 56 /*Now, one must find the eigen vectors. For that we use the following property of the inner product 57 * <Ax,y> = <x,tAy> 58 * Here, M'(M-lambda*Id) is symmetrical, which gives: 59 * ∀(x,y)∈R²xR² <M'x,y> = <M'y,x> 60 * And we have the following: 61 * if y∈Ker(M'), ∀x∈R² <M'x,y> = <x,M'y> = 0 62 * We have shown that 63 * Im(M') ⊥ Ker(M') 64 * 65 * To find the eigen vectors of M, we only have to find two vectors 66 * of the image of M' and take their perpendicular as long as they are 67 * not 0. 68 * To do that, we take the images (1,0) and (0,1): 69 * x1 = (a11 - lambda) x2 = a21 70 * y1 = a21 y2 = (a22-lambda) 71 * 72 * We take the vector that has the larger norm and take its perpendicular.*/ 73 74 double norm1 = (a11-lambda1)*(a11-lambda1) + a21*a21; 75 double norm2 = a21*a21 + (a22-lambda1)*(a22-lambda1); 76 77 if (norm2<norm1){ 78 norm1=sqrt(norm1); 79 v.x = - a21/norm1; 80 v.y = (a11-lambda1)/norm1; 81 } 82 else{ 83 norm2=sqrt(norm2); 84 v.x = - (a22-lambda1)/norm2; 85 v.y = a21/norm2; 86 } 28 87 } 29 88 30 //general case: construction of 2 eigen vectors31 else {32 delta = sqrt(delta);33 lambda1 = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0;34 double v0 = a11-lambda1, v1 = a21,v2 = a22 - lambda1;35 double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;36 if(s1 < s0)37 s0=sqrt(s0),v.x=v1/s0,v.y=-v0/s0;38 else39 s1=sqrt(s1),v.x=v2/s1,v.y=-v1/s1;40 };41 89 } 42 90 /*}}}1*/ -
issm/trunk/src/c/objects/Bamg/Mesh.cpp
r5340 r5397 21 21 if(bamggeom->Edges) { 22 22 Gh.ReadGeometry(bamggeom,bamgopts); 23 Gh. AfterRead();23 Gh.PostRead(); 24 24 } 25 25 … … 32 32 printf("WARNING: mesh present but no geometry found. Reconstructing...\n"); 33 33 BuildGeometryFromMesh(bamgopts); 34 Gh. AfterRead();34 Gh.PostRead(); 35 35 } 36 36 … … 141 141 //double cutoffradian = 10.0/180.0*Pi; 142 142 BuildGeometryFromMesh(bamgopts); 143 Gh. AfterRead();143 Gh.PostRead(); 144 144 SetIntCoor(); 145 145 ReconstructExistingMesh(); … … 300 300 BuildGeometryFromMesh(); 301 301 if (verbose) printf("Completing geometry\n"); 302 Gh. AfterRead();302 Gh.PostRead(); 303 303 } 304 304 /*}}}1*/ … … 5319 5319 int jedge=bcurve[icurve]%2; 5320 5320 5321 /*Skip if we are on a equi curve (duplicate)*/5322 if(!Gh.curves[icurve].master) continue;5323 5324 5321 /*Get edge of Bth with index iedge*/ 5325 5322 Edge &ei = BTh.edges[iedge]; -
issm/trunk/src/c/objects/Bamg/Metric.cpp
r5151 r5397 269 269 for(i=0;i<2;i++){ 270 270 /*Now, one must find the eigen vectors. For that we use the 271 * following property of the scalareproduct271 * following property of the inner product 272 272 * (Ax,b) = transp(b) Ax = transp(x) transp(A) b 273 273 * = (transp(A) b ,x)
Note:
See TracChangeset
for help on using the changeset viewer.
