Changeset 5605


Ignore:
Timestamp:
08/27/10 08:37:12 (15 years ago)
Author:
Mathieu Morlighem
Message:

As usual

Location:
issm/trunk/src/c/objects/Bamg
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Bamg/Mesh.cpp

    r5584 r5605  
    10591059        void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3) {
    10601060                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)*/
    1061                 // -------------------------------------------
     1061                // -------------------------------
    10621062                //             s2
    1063                 //                                            !
    1064                 //             /|\                            !
    1065                 //            / | \                           !
    1066                 //           /  |  \                          !
    1067                 //    tt1   /   |   \ tt0                     !
    1068                 //         /    |s   \                        !
    1069                 //        /     .     \                       !
    1070                 //       /  .      `   \                      !
    1071                 //      / .           ` \                     !
    1072                 //      ----------------                      !
     1063                //                               !
     1064                //             /|\               !
     1065                //            / | \              !
     1066                //           /  |  \             !
     1067                //    tt1   /   |   \ tt0        !
     1068                //         /    |s   \           !
     1069                //        /     .     \          !
     1070                //       /  .      `   \         !
     1071                //      / .           ` \        !
     1072                //      ----------------         !
    10731073                //   s0       tt2       s1
    1074                 //--------------------------------------------
    1075 
    1076                 //the three triangles tt
    1077                 Triangle* tt[3];
    1078                 //three vertices of t
    1079                 BamgVertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2];
    1080                 //three determinants
    1081                 Icoor2 det3local[3];
    1082                 // number of zero det3
    1083                 register int nbzerodet =0;
    1084                 // izerodet = egde contening the vertex s
    1085                 register int izerodet=-1,iedge;
     1074                //-------------------------------
     1075
     1076                /*Intermediaries*/
     1077                Triangle* tt[3];       //the three triangles
     1078                Icoor2 det3local[3];   //three determinants (integer)
     1079                int nbzerodet =0;      //number of zeros in det3
     1080                int izerodet=-1;       //egde containing the vertex s
     1081                int iedge;
     1082
     1083                /*three vertices of t*/
     1084                BamgVertex &s0=(*t)[0];
     1085                BamgVertex &s1=(*t)[1];
     1086                BamgVertex &s2=(*t)[2];
     1087
    10861088                //determinant of t
    1087                 Icoor2 detOld = t->det;
    1088 
    1089                 // infinitevertexpos = order of the infinite vertex (NULL)
    1090                 // if no infinite vertex (NULL) infinitevertexpos=-1
    1091                 // else if v_i is infinite, infinitevertexpos=i
    1092                 int infinitevertexpos = &s0 ?  ((  &s1 ? ( &s2  ? -1 : 2) : 1  )) : 0;
     1089                Icoor2 detOld=t->det;
     1090
     1091                /* infvertexindex = index of the infinite vertex (NULL)
     1092                        if no infinite vertex (NULL) infvertexindex=-1
     1093                        else if v_i is infinite, infvertexindex=i*/
     1094                int infvertexindex = &s0 ?  ((  &s1 ? ( &s2  ? -1 : 2) : 1  )) : 0;
    10931095
    10941096                //some checks
    1095                 if (( infinitevertexpos <0 ) && (detOld <0) ||  ( infinitevertexpos >=0  ) && (detOld >0) ){
    1096                         ISSMERROR("bug in Mesh::Add, bad configuration");
     1097                if (( infvertexindex <0 ) && (detOld <0) ||  ( infvertexindex >=0  ) && (detOld >0) ){
     1098                        ISSMERROR("inconsistent configuration (Contact ISSM developers)");
    10971099                }
    10981100
    10991101                // if det3 does not exist, build it
    1100                 if (!det3) {
     1102                if (!det3){
    11011103                        //allocate
    11021104                        det3 = det3local;
    11031105                        //if no infinite vertex
    1104                         if (infinitevertexpos<0 ) {
     1106                        if (infvertexindex<0 ) {
    11051107                                det3[0]=bamg::det(s ,s1,s2);
    11061108                                det3[1]=bamg::det(s0,s ,s2);
    11071109                                det3[2]=bamg::det(s0,s1,s );}
    11081110                        else {
    1109                                 // one of &s1  &s2  &s0 is NULL so (&si || &sj) <=> !&sk
     1111                                // one of &s1  &s2  &s0 is NULL
    11101112                                det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ;
    11111113                                det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ;
    … …  
    11191121
    11201122                //if nbzerodet>0, point s is on an egde or on a vertex
    1121                 if  (nbzerodet >0 ){
    1122                         if (nbzerodet == 1) {
     1123                if  (nbzerodet>0){
     1124                        /*s is on an edge*/
     1125                        if (nbzerodet==1) {
    11231126                                iedge = OppositeEdge[izerodet];
    11241127                                AdjacentTriangle ta = t->Adj(iedge);
    11251128
    1126                                 // the point is on the edge
    1127                                 // if the point is one the boundary
    1128                                 // add the point in outside part
    1129                                 if ( t->det >=0) { // inside triangle
    1130                                         if ((( Triangle *) ta)->det < 0 ) {
     1129                                /*if the point is one the boundary
     1130                                  add the point in outside part */
     1131                                if (t->det>=0){ // inside triangle
     1132                                        if (((Triangle*)ta)->det<0 ) {
    11311133                                                // add in outside triangle
    1132                                                 AddVertex(s,( Triangle *) ta);
     1134                                                AddVertex(s,( Triangle *)ta);
    11331135                                                return;
    11341136                                        }
    11351137                                }
    11361138                        }
    1137                         else {
    1138                                 //t->Echo();
    1139                                 printf("\nproblem while trying to add:\n");
    1140                                 s.Echo();
    1141                                 ISSMERROR("Bug in Mesh::Add points duplicated %i times",nbzerodet);
     1139                        else{
     1140                                ISSMERROR("Cannot add a vertex more than once. Check duplicates");
    11421141                        }
    11431142                }
    11441143
    11451144                // remove de MarkUnSwap edge
    1146                 t->SetUnMarkUnSwap(0);     
    1147                 t->SetUnMarkUnSwap(1);     
     1145                t->SetUnMarkUnSwap(0);
     1146                t->SetUnMarkUnSwap(1);
    11481147                t->SetUnMarkUnSwap(2);
    11491148
    … …  
    27902789                 * the initial simple mesh from this edge and 2 boundary triangles*/
    27912790
    2792                 BamgVertex *  v0=orderedvertices[0], *v1=orderedvertices[1];
     2791                BamgVertex *v0=orderedvertices[0], *v1=orderedvertices[1];
    27932792
    27942793                nbt = 2;
    2795                 triangles[0](0) = NULL; //infinite vertex
     2794                triangles[0](0) = NULL;//infinite vertex
    27962795                triangles[0](1) = v0;
    27972796                triangles[0](2) = v1;
    … …  
    28232822
    28242823                /*Now, add the vertices One by One*/
    2825 
    28262824                long NbSwap=0;
    28272825                if (verbose>3) printf("   Begining of insertion process...\n");
    … …  
    28332831
    28342832                        //Find the triangle in which newvertex is located
    2835                         Icoor2 dete[3];
    2836                         Triangle* tcvi = TriangleFindFromCoord(newvertex->i,dete); //(newvertex->i = integer coordinates)
     2833                        Icoor2 det3[3];
     2834                        Triangle* tcvi = TriangleFindFromCoord(newvertex->i,det3); //(newvertex->i = integer coordinates)
    28372835
    28382836                        //Add newvertex to the quadtree
    … …  
    28402838
    28412839                        //Add newvertex to the existing mesh
    2842                         AddVertex(*newvertex,tcvi,dete);
     2840                        AddVertex(*newvertex,tcvi,det3);
    28432841
    28442842                        //Make the mesh Delaunay around newvertex by swaping the edges
    … …  
    28812879                long i;
    28822880                long NbSwap=0;
    2883                 Icoor2 dete[3];
     2881                Icoor2 det3[3];
    28842882
    28852883                //number of new points
    … …  
    29252923                                }
    29262924                                vj.ReferenceNumber=0;
    2927                                 Triangle *tcvj=TriangleFindFromCoord(vj.i,dete);
     2925                                Triangle *tcvj=TriangleFindFromCoord(vj.i,det3);
    29282926                                if (tcvj && !tcvj->link){
    29292927                                        tcvj->Echo();
    … …  
    29312929                                }
    29322930                                quadtree->Add(vj);
    2933                                 AddVertex(vj,tcvj,dete);
     2931                                AddVertex(vj,tcvj,det3);
    29342932                                NbSwap += vj.Optim(1);         
    29352933                                iv++;
    … …  
    35693567        for (int icount=2; icount<nbvb; icount++) {
    35703568                BamgVertex *vi  = orderedvertices[icount];
    3571                 Icoor2 dete[3];
    3572                 Triangle *tcvi = TriangleFindFromCoord(vi->i,dete);
     3569                Icoor2 det3[3];
     3570                Triangle *tcvi = TriangleFindFromCoord(vi->i,det3);
    35733571                quadtree->Add(*vi);
    3574                 AddVertex(*vi,tcvi,dete);
     3572                AddVertex(*vi,tcvi,det3);
    35753573                NbSwap += vi->Optim(1,1);
    35763574        }
    … …  
    46824680                long  iv = nbvold;
    46834681                long NbSwap = 0;
    4684                 Icoor2 dete[3]; 
     4682                Icoor2 det3[3]; 
    46854683                for (int i=nbvold;i<nbv;i++) {// for all the new point
    46864684                        BamgVertex & vi = vertices[i];
    … …  
    46914689                        vi.ReferenceNumber=0;
    46924690                        vi.DirOfSearch =NoDirOfSearch;
    4693                         Triangle *tcvi = TriangleFindFromCoord(vi.i,dete);
     4691                        Triangle *tcvi = TriangleFindFromCoord(vi.i,det3);
    46944692                        if (tcvi && !tcvi->link) {
    46954693                                printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
    … …  
    47004698                                ISSMERROR("!tcvi || tcvi->det < 0");
    47014699                        }
    4702                         AddVertex(vi,tcvi,dete);
     4700                        AddVertex(vi,tcvi,det3);
    47034701                        NbSwap += vi.Optim(1);         
    47044702                        iv++;
    … …  
    47284726/*}}}1*/
    47294727/*FUNCTION Mesh::TriangleFindFromCoord{{{1*/
    4730 Triangle * Mesh::TriangleFindFromCoord(const I2 & B,Icoor2 dete[3], Triangle *tstart) const {
     4728Triangle * Mesh::TriangleFindFromCoord(const I2 & B,Icoor2 det3[3], Triangle *tstart) const {
    47314729        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/
    47324730
    … …  
    47694767                ISSMASSERT(k0>=0);// k0 the NULL vertex
    47704768                int k1=NextVertex[k0],k2=PreviousVertex[k0];
    4771                 dete[k0]=det(B,(*t)[k1],(*t)[k2]);
    4772                 dete[k1]=dete[k2]=-1;     
    4773                 if (dete[k0] > 0) // outside B
     4769                det3[k0]=det(B,(*t)[k1],(*t)[k2]);
     4770                det3[k1]=det3[k2]=-1;     
     4771                if (det3[k0] > 0) // outside B
    47744772                 return t;
    47754773                t = t->TriangleAdj(OppositeEdge[k0]);
    … …  
    47874785
    47884786                j= OppositeVertex[jj];
    4789                 dete[j] = detop;  //det(*b,*s1,*s2);
     4787                det3[j] = detop;  //det(*b,*s1,*s2);
    47904788                jn = NextVertex[j];
    47914789                jp = PreviousVertex[j];
    4792                 dete[jp]= det(*(*t)(j),*(*t)(jn),B);
    4793                 dete[jn] = t->det-dete[j] -dete[jp];
    4794 
    4795                 // count the number k of  dete <0
     4790                det3[jp]= det(*(*t)(j),*(*t)(jn),B);
     4791                det3[jn] = t->det-det3[j] -det3[jp];
     4792
     4793                // count the number k of  det3 <0
    47964794                int k=0,ii[3];
    4797                 if (dete[0] < 0 ) ii[k++]=0;
    4798                 if (dete[1] < 0 ) ii[k++]=1;
    4799                 if (dete[2] < 0 ) ii[k++]=2;
     4795                if (det3[0] < 0 ) ii[k++]=0;
     4796                if (det3[1] < 0 ) ii[k++]=1;
     4797                if (det3[2] < 0 ) ii[k++]=2;
    48004798                // 0 => ok
    48014799                // 1 => go in way 1
    … …  
    48034801
    48044802                if (k==0) break;
    4805                 if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]);
    4806                 if ( k>=3){
    4807                         ISSMERROR("k>=3");
    4808                 }
     4803                if (k==2 && BinaryRand()) Exchange(ii[0],ii[1]);
     4804                ISSMASSERT(k<3);
    48094805                AdjacentTriangle t1 = t->Adj(jj=ii[0]);
    48104806                if ((t1.det() < 0 ) && (k == 2))
    … …  
    48124808                t=t1;
    48134809                j=t1;// for optimisation we now the -det[OppositeVertex[j]];
    4814                 detop = -dete[OppositeVertex[jj]];
     4810                detop = -det3[OppositeVertex[jj]];
    48154811                jj = j;
    48164812        }
    48174813
    48184814        if (t->det<0) // outside triangle
    4819          dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;
     4815         det3[0]=det3[1]=det3[2]=-1,det3[OppositeVertex[jj]]=detop;
    48204816        return t;
    48214817}

  • issm/trunk/src/c/objects/Bamg/Triangle.cpp

    r5460 r5605  
    259259        /*FUNCTION Triangle::SetAdj2{{{1*/
    260260        void Triangle::SetAdj2(short a,Triangle *t,short aat){
    261                 adj[a]=t;    //the adjacent triangle to the edge a is t
    262                 AdjEdgeIndex[a]=aat; //position of the edge in the adjacent triangle
    263                 if(t) { //if t!=NULL add adjacent triangle to t (this)
     261                /*For current triangle:
     262                 * - a is the index of the edge were the adjency is set (in [0 2])
     263                 * - t is the adjacent triangle
     264                 * - aat is the index of the same edge in the adjacent triangle*/
     265                adj[a]=t;
     266                AdjEdgeIndex[a]=aat;
     267                if(t){ //if t!=NULL add adjacent triangle to t (this)
    264268                        t->adj[aat]=this;
    265269                        t->AdjEdgeIndex[aat]=a;
Note: See TracChangeset for help on using the changeset viewer.