Changeset 5339

Timestamp:

08/18/10 08:32:08 (15 years ago)

Author:

Mathieu Morlighem

Message:

As usual

File:

• issm/trunk/src/c/objects/Bamg/Geometry.cpp (modified) (9 diffs)

issm/trunk/src/c/objects/Bamg/Geometry.cpp

@@ -478 +478 @@
                  * a vertex v?
                  * To do so, we use the same chaining algorithm but there is a difficulty
-                 * coming from the fact that F we have a couple of vertices and not one 
-                 * vertices.
-                 * To overcome this difficulty, we use a global indices exactly like in 
+                 * coming from the fact that for F we have a couple of vertices and not one 
+                 * vertex.
+                 * To overcome this difficulty, we use a global indexing exactly like in 
                  * C/C++ so that
                  * a member of a 2-column-table can be described by one index p=i*2+j
                  * i=(int)p/2 line number of p
-                 * j=p%2       column number of p
+                 * j=p%2      column number of p
                  *
                  * Algorithm:
@@ -506 +506 @@
                         double lv10=Norme2(v10);
                         //check that its length is not 0
-                        if(lv10==0) {
-                                ISSMERROR("Length of edge %i is 0",i);
-                        }
+                        if(lv10==0)ISSMERROR("Length of edge %i is 0",i);
                         //compute angle in [-Pi Pi]
                         eangle[i] = atan2(v10.y,v10.x);
@@ -521 +519 @@
                 //sort head_v by order of increasing edges angle
                 for (i=0;i<nbv;i++) {
-                        int exch=1, ord=0;      
+                        int exch=1,ord=0;      
 
                         //exchange vertices position in head_v and next_p till tey are sorted
                         while (exch){
-                                long *p=head_v+i;               // pointer toward head_v[vertex i]
-                                long *po=p;                     // copy of pointer p
-                                long  n=*p;                     // next value of edge holding i
+                                long *p=head_v+i;               
+                                long *po=p;                     
+                                long  n=*p;                     
                                 register float angleold=-1000 ; // angle = - infinity
                                 ord=0; exch=0;
 
-                                // loop over the edges that contain the vertex i
+                                // loop over the edges that contain the vertex i (till -1)
                                 while (n >=0){
                                         ord++;
-                                        register long  i1=n/2;       // i1 = floor (n/2)
-                                        register long  j1=n%2;       // j1 = 1 if n is odd
-                                        register long* pn=next_p+n;  // pointer to next_p[n]
-
-                                        //  n = next_p[n] = position in edge of next vertex i
+                                        long  i1=n/2;       // i1 = floor (n/2)      -> Edge number
+                                        long  j1=n%2;       // j1 = 1 if n is odd    -> Vertex index for this edge (0 or 1)
+                                        long* pn=next_p+n;
+
+                                        //Next vertex index
                                         n=*pn;                       
 
@@ -562 +560 @@
                         }
 
-                        // angular test on current vertex to guess whether it is a corner (ord = number of edges horlding i)
-                        if(ord == 2) { 
+                        // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i)
+                        if(ord==2) { 
                                 long  n1 = head_v[i];
                                 long  n2 = next_p[n1];
@@ -574 +572 @@
                                         vertices[i].SetCorner() ; 
                                 }
-                                // if the ReferenceNumber a changing then is     SetRequired();
+                                // if the edge type/referencenumber a changing then is SetRequired();
                                 if (edges[i1].type != edges[i2].type || edges[i1].Required()){
                                         vertices[i].SetRequired();
@@ -586 +584 @@
                         }
 
-                        // close the list around the vertex 
+                        /*close the list around the vertex to have a circular loop*/
                         long no=-1, ne = head_v[i];
                         while (ne >=0) ne = next_p[no=ne];        
                         if(no>=0) next_p[no] = head_v[i];
-                        // now the list around the vertex is circular
-                }
-
+                }
+
+                /*Check that the list makes sense (we have all the time the same vertex)
+                 * and build adjacent edges*/
                 k =0;
                 for (i=0;i<nbe;i++){
                         for (j=0;j<2;j++){
+
                                 long n1 = next_p[k++]; 
                                 long i1 = n1/2 ,j1=n1%2;
-                                if( edges[i1].v[j1] != edges[i].v[j]) {
-                                        ISSMERROR("Bug Adj edge");
-                                }
+
+                                if( edges[i1].v[j1] != edges[i].v[j]) ISSMERROR("Problem while processing edges: check the edge list");
+
                                 edges[i1].Adj[j1] = edges + i;
                                 edges[i1].DirAdj[j1] = j;
@@ -608 +608 @@
                 // generation of  all the tangents
                 for (i=0;i<nbe;i++) {
-                        //Get AB = vertex1->vertex2
-                        R2    AB =edges[i].v[1]->r -edges[i].v[0]->r;        
-                        //Get length of AB
-                        double lAB=Norme2(AB);
-                        //initialize tangent
-                        double ltg2[2]={0.0};
+                        R2    AB =edges[i].v[1]->r -edges[i].v[0]->r;// AB = vertex0 -> vertex1
+                        double lAB=Norme2(AB);                       // Get length of AB
+                        double ltg2[2]={0.0};                        // initialize tangent
 
                         //loop over the 2 vertices of the edge
-                        for (jj=0;jj<2;jj++) {
-                                R2     tg =edges[i].tg[jj];
+                        for (j=0;j<2;j++) {
+                                R2     tg =edges[i].tg[j];
                                 double ltg=Norme2(tg);
 
                                 //by default, tangent=[0 0]
-                                if(ltg == 0) {
+                                if(ltg==0){
                                         //if the current vertex of the edge is not a corner
-                                        if(!edges[i].v[jj]->Corner()){
-                                                tg = edges[i].v[1-jj]->r - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r; //previous point and next point on curve
+                                        if(!edges[i].v[j]->Corner()){
+                                                /*The tangent is set as the vector between the
+                                                 * previous and next vertices connected to current vertex
+                                                 * normed by the edge length*/
+                                                tg = edges[i].v[1-j]->r - edges[i].Adj[j]->v[1-edges[i].DirAdj[j]]->r;
                                                 ltg= Norme2(tg);
                                                 tg = tg *(lAB/ltg);
                                                 ltg= lAB;
                                         }
-                                        //else:  a Corner with no tangent => nothing to do    
+                                        //else:  a Corner no tangent => nothing to do    
                                 }
                                 else{
@@ -636 +636 @@
                                 }
 
-                                ltg2[jj] = ltg;
-
-                                if ((tg,AB)<0){
-                                        tg = -tg;
-                                }
-
-                                edges[i].tg[jj] = tg;
+                                ltg2[j] = ltg;
+
+                                if ((tg,AB)<0) tg = -tg;
+
+                                edges[i].tg[j]=tg;
                         }
                         if (ltg2[0]!=0) edges[i].SetTgA();
@@ -649 +647 @@
 
                 for (int step=0;step<2;step++){
+
                         for (i=0;i<nbe;i++) edges[i].SetUnMark();
+
                         nbcurves = 0;
                         long  nbgem=0;
-                        for (int level=0;level < 2 && nbgem != nbe;level++)
-                         for (i=0;i<nbe;i++) {
-                                 GeometricalEdge & ei = edges[i];   
-                                 for(jj=0;jj<2;jj++){
-                                         if (!ei.Mark() && (level || ei[jj].Required())) { 
-                                                 // warning ei.Mark() can be change in loop for(jj=0;jj<2;jj++) 
-                                                 int k0=jj,k1;
-                                                 GeometricalEdge *e = & ei;
-                                                 GeometricalVertex *a=(*e)(k0); // begin 
-                                                 if(curves){
-                                                         curves[nbcurves].be=e;
-                                                         curves[nbcurves].kb=k0;
-                                                 }
-                                                 int nee=0;
-                                                 for(;;) { 
-                                                         nee++;
-                                                         k1 = 1-k0; // next vertex of the edge 
-                                                         e->SetMark();
-                                                         nbgem++;
-                                                         e->CurveNumber=nbcurves;
-                                                         if(curves) {
-                                                                 curves[nbcurves].ee=e;
-                                                                 curves[nbcurves].ke=k1;
-                                                         }
-
-                                                         GeometricalVertex *b=(*e)(k1);
-                                                         if (a == b ||  b->Required() ) break;
-                                                         k0 = e->DirAdj[k1];//  vertex in next edge
-                                                         e = e->Adj[k1]; // next edge
-                                                 }
-                                                 nbcurves++;
-                                                 if(level) a->SetRequired();
-                                         }
-                                 }
-                         } 
+
+                        for (int level=0;level < 2 && nbgem != nbe;level++){
+                                for (i=0;i<nbe;i++){
+                                        GeometricalEdge & ei = edges[i];   
+                                        for(j=0;j<2;j++){
+                                                if (!ei.Mark() && (level || ei[j].Required())) { 
+                                                        // warning ei.Mark() can be change in loop for(j=0;j<2;j++) 
+                                                        int k0=j,k1;
+                                                        GeometricalEdge *e = & ei;
+                                                        GeometricalVertex *a=(*e)(k0); // begin 
+                                                        if(curves){
+                                                                curves[nbcurves].be=e;
+                                                                curves[nbcurves].kb=k0;
+                                                        }
+                                                        int nee=0;
+                                                        for(;;) { 
+                                                                nee++;
+                                                                k1 = 1-k0; // next vertex of the edge 
+                                                                e->SetMark();
+                                                                nbgem++;
+                                                                e->CurveNumber=nbcurves;
+                                                                if(curves) {
+                                                                        curves[nbcurves].ee=e;
+                                                                        curves[nbcurves].ke=k1;
+                                                                }
+
+                                                                GeometricalVertex *b=(*e)(k1);
+                                                                if (a == b ||  b->Required() ) break;
+                                                                k0 = e->DirAdj[k1];//  vertex in next edge
+                                                                e = e->Adj[k1]; // next edge
+                                                        }
+                                                        nbcurves++;
+                                                        if(level) a->SetRequired();
+                                                }
+                                        }
+                                } 
+                        }
                         ISSMASSERT(nbgem && nbe);
                         if(step==0){