Ignore:
Timestamp:
10/28/13 14:43:03 (11 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 16554

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c/bamg/Mesh.cpp

    r15219 r16560  
    122122                  }
    123123                  for (i=0;i<Tho.nbt;i++)
    124                         if(  reft[i] >=0 && flag[i])
    125                           {
     124                        if(reft[i] >=0 && flag[i]){
    126125                                const Triangle & t = Tho.triangles[i];
    127126                                int i0 = Tho.GetId(t[0]);
    128127                                int i1 = Tho.GetId(t[1]);
    129128                                int i2 = Tho.GetId(t[2]);
    130                                 if (i0<0 || i1<0 || i2<0){
     129                                if(i0<0 || i1<0 || i2<0){
    131130                                        delete [] refv;
    132131                                        _error_("i0<0 || i1<0 || i2< 0");
    133132                                }
    134                                 if (i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv){
     133                                if(i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv){
     134                                        delete [] refv;
    135135                                        _error_("i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv");
    136136                                }
    … …  
    139139                                nbt++;           
    140140                          }
    141                   if (kt!=nbt){
    142                           _error_("kt!=nbt");
    143                   }
    144141                  if (nbt==0 && nbv==0) {
    145142                          _error_("All triangles have been removed");
    … …  
    148145                  delete [] reft;
    149146                  delete [] refv;
    150                   //double cutoffradian = 10.0/180.0*Pi;
    151147                  BuildGeometryFromMesh(bamgopts);
    152148                  Gh.PostRead();
    … …  
    154150                  ReconstructExistingMesh();
    155151
    156                   if (!nbsubdomains){
    157                           _error_("nbsubdomains==0");
    158                   }
    159                   if (!subdomains[0].head || !subdomains[0].head->link){
    160                           _error_("!subdomains[0].head || !subdomains[0].head->link");
    161                   }
    162 
     152                  /*Final checks*/
     153                  _assert_(kt==nbt);
     154                  _assert_(nbsubdomains);
     155                  _assert_(subdomains[0].head && subdomains[0].head->link);
    163156          }
    164157        /*}}}*/
    … …  
    269262        void Mesh::ReadMesh(int* index,double* x,double* y,int nods,int nels){
    270263
    271                 double Hmin = HUGE_VAL;// the infinie value
    272264                long i1,i2,i3;
    273265                long i;
    … …  
    523515                                head=(int)bamgmesh->SubDomains[i*3+1]-1;//C indexing
    524516                                direction=(int)bamgmesh->SubDomains[i*3+2];
    525                                 if (i3!=23) _error_("Bad Subdomain definition: first number should be 3");
     517                                if (i3!=3) _error_("Bad Subdomain definition: first number should be 3");
    526518                                if (head<0 || head>=nbt) _error_("Bad Subdomain definition: head should in [1 " << nbt << "] (triangle number)");
    527519                                subdomains[i].head = triangles+head;
     520                                subdomains[i].direction = direction;
     521                                subdomains[i].ReferenceNumber = i3;
    528522                        }
    529523                }
    … …  
    10301024                                Triangle &t=triangles[i];
    10311025                                if (t.det>0 && !(t.Hidden(0)||t.Hidden(1) || t.Hidden(2) )){
     1026                                        /*Remove triangles that have a bad aspect ratio*/
    10321027                                        //if(t.Anisotropy()<2 & t.Length()<1.e+5){
    10331028                                                index[num*3+0]=GetId(t[0])+1; //back to M indexing
    … …  
    10521047
    10531048                /*Get options*/
    1054                 int    verbose=bamgopts->verbose;
    1055                 double anisomax =bamgopts->anisomax;
    1056                 double errg     =bamgopts->errg;
     1049                double anisomax = bamgopts->anisomax;
     1050                double errg     = bamgopts->errg;
    10571051
    10581052                double ss[2]={0.00001,0.99999};
    … …  
    10621056
    10631057                //check that hmax is positive
    1064                 if (hmax<=0){
    1065                         _error_("hmax<=0");
    1066                 }
     1058                _assert_(hmax>0);
    10671059
    10681060                //errC cannot be higher than 1
    1069                 if (errC>1) errC=1;
     1061                if(errC>1) errC=1;
    10701062
    10711063                //Set all vertices to "on"
    … …  
    10731065
    10741066                //loop over all the vertices on edges
    1075                 for (int  i=0;i<nbe;i++){
     1067                for(int  i=0;i<nbe;i++){
    10761068                        for (int j=0;j<2;j++){
    10771069
    … …  
    11601152
    11611153                //some checks
    1162                 if (( infvertexindex <0 ) && (detOld <0) ||  ( infvertexindex >=0  ) && (detOld >0) ){
     1154                if(((infvertexindex <0 ) && (detOld <0)) ||  ((infvertexindex >=0) && (detOld >0)) ){
    11631155                        _error_("inconsistent configuration (Contact ISSM developers)");
    11641156                }
    … …  
    17271719                        long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
    17281720                        k = edge4->SortAndFind(i0,i1);
    1729                         if(k>=0){
    1730                                 subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
    1731                                 subdomains[i].edge = edges+k;
    1732                                 Gh.subdomains[i].edge = Gh.edges + k;
    1733                                 Gh.subdomains[i].direction  =  subdomains[i].direction;
    1734                                 Gh.subdomains[i].ReferenceNumber =  subdomains[i].ReferenceNumber;
    1735                         }
    1736                         else
    1737                          _error_("%i should be >=0");
    1738                   }
     1721                        _assert_(k>=0);
     1722                        subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
     1723                        subdomains[i].edge = edges+k;
     1724                        Gh.subdomains[i].edge = Gh.edges + k;
     1725                        Gh.subdomains[i].direction  =  subdomains[i].direction;
     1726                        Gh.subdomains[i].ReferenceNumber =  subdomains[i].ReferenceNumber;
     1727                }
    17391728
    17401729                delete edge4;
    … …  
    21872176                                // correction of second derivative
    21882177                                // by a laplacien
    2189                                 double* d2[3] = {dxdx, dxdy, dydy};
    21902178                                double* dd;
    21912179                                for (int xy = 0;xy<3;xy++) {
    2192                                         dd = d2[xy];
     2180                                        if      (xy==0) dd=dxdx;
     2181                                        else if (xy==1) dd=dxdy;
     2182                                        else if (xy==2) dd=dydy;
     2183                                        else    _error_("not supported yet");
    21932184                                        // do leat 2 iteration for boundary problem
    21942185                                        for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
    … …  
    26202611
    26212612                          }
    2622                         else
    2623                           { // find the head for all sub domaine
     2613                        else{ // find the head for all subdomains
    26242614                                if (Gh.nbsubdomains != nbsubdomains && subdomains)
    26252615                                 delete [] subdomains, subdomains=0;
    … …  
    26272617                                 subdomains = new SubDomain[ Gh.nbsubdomains];
    26282618                                nbsubdomains =Gh.nbsubdomains;
    2629                                 long err=0;
    26302619                                CreateSingleVertexToTriangleConnectivity();
    26312620                                long * mark = new long[nbt];
    … …  
    26382627                                        GeomEdge &eg = *Gh.subdomains[i].edge;
    26392628                                        subdomains[i].ReferenceNumber = Gh.subdomains[i].ReferenceNumber;
    2640                                         int ssdlab = subdomains[i].ReferenceNumber;
    26412629                                        // by carefull is not easy to find a edge create from a GeomEdge
    26422630                                        // see routine MakeGeomEdgeToEdge
    … …  
    27732761                        vertices=new BamgVertex[maxnbv];
    27742762                        _assert_(vertices);
    2775                         orderedvertices=new (BamgVertex* [maxnbv]);
     2763                        orderedvertices=new BamgVertex* [maxnbv];
    27762764                        _assert_(orderedvertices);
    27772765                        triangles=new Triangle[maxnbt];
    … …  
    30153003        }
    30163004        /*}}}*/
    3017         /*FUNCTION Mesh::isCracked{{{*/
    3018         int Mesh::isCracked() const {
    3019                 return NbCrackedVertices != 0;
    3020         }
    3021         /*}}}*/
    30223005        /*FUNCTION Mesh::MakeGeomEdgeToEdge{{{*/
    30233006        Edge** Mesh::MakeGeomEdgeToEdge() {
    … …  
    30273010                        _error_("!Gh.nbe");
    30283011                }
    3029                 Edge **e= new (Edge* [Gh.nbe]);
     3012                Edge **e= new Edge* [Gh.nbe];
    30303013
    30313014                long i;
    … …  
    31083091        void Mesh::MakeBamgQuadtree() { 
    31093092                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)*/
    3110 
    3111                 long int verbose=0;
    3112                 if (  !quadtree )  quadtree = new BamgQuadtree(this);
    3113 
     3093                if(!quadtree) quadtree = new BamgQuadtree(this);
    31143094        }
    31153095        /*}}}*/
    … …  
    31353115                        for (int j=0;j<3;j++){
    31363116                                Triangle &tt = *t.TriangleAdj(j);
    3137                                 if ( ! &tt ||  it < GetId(tt) && ( tt.link || t.link)){
     3117                                if ( (!&tt ||  it < GetId(tt)) && ( tt.link || t.link)){
    31383118                                        BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
    31393119                                        BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
    … …  
    33133293                        }
    33143294
    3315                 } while (nbv!=nbvold);
    3316 
    3317                 delete []  first_np_or_next_t;
    3318 
    3319                 long NbSwapf =0,NbSwp;
    3320 
    3321                 NbSwp = NbSwapf;
    3322                 for (i=0;i<nbv;i++)
    3323                  NbSwapf += vertices[i].Optim(0);
    3324                 NbTSwap +=  NbSwapf ;
    3325         }
    3326         /*}}}*/
     3295                }while(nbv!=nbvold);
     3296                delete [] first_np_or_next_t;
     3297
     3298                long NbSwapf =0;
     3299                for(i=0;i<nbv;i++) NbSwapf += vertices[i].Optim(0);
     3300        }/*}}}*/
    33273301        /*FUNCTION Mesh::ProjectOnCurve{{{*/
    33283302        GeomEdge*   Mesh::ProjectOnCurve( Edge & BhAB, BamgVertex &  vA, BamgVertex & vB,
    … …  
    37653739                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
    37663740
    3767                 long int verbose=0;
    37683741                long *renu= new long[nbt];
    37693742                register Triangle *t0,*t,*te=triangles+nbt;
    … …  
    38303803        }
    38313804        /*}}}*/
    3832         /*FUNCTION Mesh::VerticesRenumber{{{*/
    3833         void Mesh::VerticesRenumber(long * renu) {
    3834                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
    3835 
    3836                 // warning be carfull because pointer
    3837                 // from on mesh to over mesh
    3838                 //  --  so do ReNumbering at the beginning
    3839                 BamgVertex * ve = vertices+nbv;
    3840                 long it,ie,i;
    3841 
    3842                 _printf_("renumbering triangles\n");
    3843                 for ( it=0;it<nbt;it++)
    3844                  triangles[it].Renumbering(vertices,ve,renu);
    3845 
    3846                 _printf_("renumbering edges\n");
    3847                 for ( ie=0;ie<nbe;ie++)
    3848                  edges[ie].Renumbering(vertices,ve,renu);
    3849 
    3850                 _printf_("renumbering vertices on geom\n");
    3851                 for (i=0;i< NbVerticesOnGeomVertex;i++)
    3852                   {
    3853                         BamgVertex *v = VerticesOnGeomVertex[i].meshvertex;
    3854                         if (v>=vertices && v < ve)
    3855                          VerticesOnGeomVertex[i].meshvertex=vertices+renu[GetId(v)];
    3856                   }
    3857 
    3858                 _printf_("renumbering vertices on edge\n");
    3859                 for (i=0;i< NbVerticesOnGeomEdge;i++)
    3860                   {
    3861                         BamgVertex *v =VerticesOnGeomEdge[i].meshvertex;
    3862                         if (v>=vertices && v < ve)
    3863                          VerticesOnGeomEdge[i].meshvertex=vertices+renu[GetId(v)];
    3864                   }
    3865 
    3866                 _printf_("renumbering vertices on Bth vertex\n");
    3867                 for (i=0;i< NbVertexOnBThVertex;i++)
    3868                   {
    3869                         BamgVertex *v=VertexOnBThVertex[i].v;
    3870                         if (v>=vertices && v < ve)
    3871                          VertexOnBThVertex[i].v=vertices+renu[GetId(v)];
    3872                   }
    3873 
    3874                 for (i=0;i< NbVertexOnBThEdge;i++)
    3875                   {
    3876                         BamgVertex *v=VertexOnBThEdge[i].v;
    3877                         if (v>=vertices && v < ve)
    3878                          VertexOnBThEdge[i].v=vertices+renu[GetId(v)];
    3879                   }
    3880 
    3881                 // move the Vertices without a copy of the array
    3882                 // be carefull not trivial code
    3883                 long j;
    3884                 for ( it=0;it<nbv;it++) // for all sub cycles of the permutation renu
    3885                  if (renu[it] >= 0) // a new sub cycle
    3886                         {
    3887                          i=it;
    3888                          BamgVertex ti=vertices[i],tj;
    3889                          while ( (j=renu[i]) >= 0){
    3890                                  // i is old, and j is new
    3891                                  renu[i] = -1-renu[i]; // mark
    3892                                  tj = vertices[j];     // save new
    3893                                  vertices[j]= ti;      // new <- old
    3894                                  i=j;     // next
    3895                                  ti = tj;
    3896                                 } 
    3897                         }
    3898                 if (quadtree){
    3899                         delete quadtree;
    3900                         quadtree = new BamgQuadtree(this);
    3901                 }
    3902 
    3903                 for ( it=0;it<nbv;it++) renu[i]= -renu[i]-1;
    3904         }
    3905         /*}}}*/
    39063805/*FUNCTION Mesh::SetIntCoor{{{*/
    39073806void Mesh::SetIntCoor(const char * strfrom) {
    … …  
    39643863
    39653864        if (number_of_errors) _error_("Fatal error: some triangles have negative areas, see above");
    3966 }
    3967 /*}}}*/
    3968 /*FUNCTION Mesh::ShowRegulaty{{{*/
    3969 void  Mesh::ShowRegulaty() const {
    3970         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr>*/
    3971 
    3972         const  double  sqrt32=sqrt(3.)*0.5;
    3973         const double  aireKh=sqrt32*0.5;
    3974         D2  Beq(1,0),Heq(0.5,sqrt32);
    3975         D2xD2 Br(D2xD2(Beq,Heq).t());
    3976         D2xD2 B1r(Br.inv());
    3977         double gammamn=1e100,hmin=1e100;
    3978         double gammamx=0,hmax=0;
    3979         double beta=1e100;
    3980         double beta0=0;
    3981         double  alpha2=0;
    3982         double area=0,Marea=0;
    3983         // double cf= double(coefIcoor);
    3984         // double cf2= 6.*cf*cf;
    3985         int nt=0;
    3986         for (int it=0;it<nbt;it++)
    3987          if ( triangles[it].link){
    3988                  Triangle &K=triangles[it];
    3989                  double  area3= Area2((R2) K[0],(R2) K[1],(R2) K[2])/6.;
    3990                  area+= area3;
    3991                  D2xD2 B_Kt(K[0],K[1],K[2]);
    3992                  D2xD2 B_K(B_Kt.t());
    3993                  D2xD2 B1K = Br*B_K.inv();
    3994                  D2xD2 BK =  B_K*B1r;
    3995                  D2xD2 B1B1 = B1K.t()*B1K;
    3996                  Metric MK(B1B1.x.x,B1B1.x.y,B1B1.y.y);
    3997                  EigenMetric VMK(MK);
    3998                  alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1));
    3999                  double betaK=0;
    4000 
    4001                  for (int j=0;j<3;j++)
    4002                         {
    4003                          double he= Norme2(R2(K[j],K[(j+1)%3]));
    4004                          hmin=Min(hmin,he);
    4005                          hmax=Max(hmax,he);
    4006                          BamgVertex & v=K[j];
    4007                          D2xD2 M((Metric)v);
    4008                          betaK += sqrt(M.det());
    4009 
    4010                          D2xD2 BMB = BK.t()*M*BK;
    4011                          Metric M1(BMB.x.x,BMB.x.y,BMB.y.y);
    4012                          EigenMetric VM1(M1);
    4013                          gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2);
    4014                          gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2);         
    4015                         }
    4016                  betaK *= area3;//  1/2 (somme sqrt(det))* area(K)
    4017                  Marea+= betaK;
    4018                  beta=min(beta,betaK);
    4019                  beta0=max(beta0,betaK);
    4020                 }   
    4021         area*=3;
    4022         gammamn=sqrt(gammamn);
    4023         gammamx=sqrt(gammamx);   
    4024         _printf_("   Adaptmesh info:\n");
    4025         _printf_("      number of triangles = " << nt << "\n");
    4026         _printf_("      hmin = " << hmin << ", hmax=" << hmax << "\n");
    4027         _printf_("      area = " << area << ", M area = " << Marea << ", M area/( |Khat| nt) = " <<  Marea/(aireKh*nt) << "\n");
    4028         _printf_("      infinite-regularity(?): min = " << gammamn << ", max = " << gammamx << "\n");
    4029         _printf_("      anisomax = " << pow(alpha2,0.5) << ", beta max = " << 1./pow(beta/aireKh,0.5) << ", min = " << 1./pow(beta0/aireKh,0.5) << "\n");
    4030 }
    4031 /*}}}*/
    4032 /*FUNCTION Mesh::ShowHistogram{{{*/
    4033 void  Mesh::ShowHistogram() const {
    4034         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ShowHistogram)*/
    4035 
    4036         const long kmax=10;
    4037         const double llmin = 0.5,llmax=2;
    4038         const double lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin);
    4039         long histo[kmax+1];
    4040         long i,it,k, nbedges =0;
    4041         for (i=0;i<=kmax;i++) histo[i]=0;
    4042         for (it=0;it<nbt;it++)
    4043          if ( triangles[it].link)
    4044                 {
    4045 
    4046                  for (int j=0;j<3;j++)
    4047                         {
    4048                          Triangle *ta = triangles[it].TriangleAdj(j);
    4049                          if ( !ta || !ta->link || GetId(ta) >= it)
    4050                                 {
    4051                                  BamgVertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]];
    4052                                  BamgVertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]];
    4053                                  if ( !& vP || !&vQ) continue;
    4054                                  R2 PQ = vQ.r - vP.r;
    4055                                  double l = log(LengthInterpole(vP,vQ,PQ));
    4056                                  nbedges++;
    4057                                  k = (int) ((l - lmin)*delta);
    4058                                  k = Min(Max(k,0L),kmax);
    4059                                  histo[k]++;
    4060                                 }
    4061                         }
    4062                 } 
    4063         _printf_(" --- Histogram of the unit mesh,  nb of edges = " << nbedges << "\n");
    4064         _printf_("      length of edge in   | %% of edge  | Nb of edges \n");
    4065         _printf_("      --------------------+-------------+-------------\n");
    4066         for   (i=0;i<=kmax;i++){
    4067                 if (i==0) _printf_( "      " << setw(10) << 0.);
    4068                 else      _printf_( "      " << setw(10) << exp(lmin+i/delta));
    4069                 if (i==kmax) _printf_("          +inf   ");
    4070                 else      _printf_( "      " << setw(10) << exp(lmin+(i+1)/delta));
    4071                 _printf_("|  " << setw(10) << (long((10000. * histo[i])/ nbedges)/100.) << " |\n");
    4072                 _printf_("  " << histo[i] << "\n");
    4073         }
    4074         _printf_("      --------------------+-------------+-------------\n");
    40753865}
    40763866/*}}}*/
    … …  
    41383928         first_np_or_next_t1[i]=-1;
    41393929        kk=0;
    4140         while (Head0>=0&& kk++<100) {
     3930        while(Head0>=0&& kk++<100){
    41413931                kch=0;
    4142                 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
     3932                for(i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
    41433933                        //  pour tous les triangles autour du sommet s
    41443934                        register Triangle* t= vertices[i].t;
    … …  
    41953985}
    41963986/*}}}*/
    4197         /*FUNCTION Mesh::SplitElement{{{*/
    4198         int  Mesh::SplitElement(int choice){
    4199                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/SplitElement)*/
    4200 
    4201                 long int verbose=0;
    4202 
    4203                 Direction NoDirOfSearch;
    4204                 const  int withBackground = &BTh != this && &BTh;
    4205 
    4206                 TrianglesRenumberBySubDomain();
    4207                 //int nswap =0;
    4208                 const long nfortria( choice ? 4 : 6);
    4209                 if(withBackground)
    4210                   {
    4211                         BTh.SetVertexFieldOn();
    4212                         SetVertexFieldOnBTh();
    4213                   }
    4214                 else
    4215                  BTh.SetVertexFieldOn();
    4216 
    4217                 long newnbt=0,newnbv=0;
    4218                 long * kedge = 0;
    4219                 long newnbq=nbq;
    4220                 long * ksplit= 0, * ksplitarray=0;
    4221                 long kkk=0;
    4222                 int ret =0;
    4223                 if (maxnbv<nbv+nbe) return 1;//   
    4224                 // 1) create  the new points by spliting the internal edges
    4225                 // set the
    4226                 long nbvold = nbv;
    4227                 long nbtold = nbt;
    4228                 long nbtoutold  = nbtout;
    4229                 long  NbEdgeOnGeom=0;
    4230                 long i;
    4231 
    4232                 nbt = nbt - nbtout; // remove all the  the ouside triangles
    4233                 long nbtsave = nbt;
    4234                 Triangle * lastT = triangles + nbt;
    4235                 for (i=0;i<nbe;i++)
    4236                  if(edges[i].GeomEdgeHook) NbEdgeOnGeom++;
    4237                 long newnbe=nbe+nbe;
    4238                 //  long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
    4239                 long newNbVerticesOnGeomEdge=NbVerticesOnGeomEdge+NbEdgeOnGeom;
    4240                 // long newNbVertexOnBThVertex=NbVertexOnBThVertex;
    4241                 long newNbVertexOnBThEdge=withBackground ? NbVertexOnBThEdge+NbEdgeOnGeom:0;
    4242 
    4243                 // do allocation for pointeur to the geometry and background
    4244                 VertexOnGeom * newVerticesOnGeomEdge = new VertexOnGeom[newNbVerticesOnGeomEdge];
    4245                 VertexOnEdge *newVertexOnBThEdge = newNbVertexOnBThEdge ?  new VertexOnEdge[newNbVertexOnBThEdge]:0;
    4246                 if (NbVerticesOnGeomEdge)
    4247                  memcpy(newVerticesOnGeomEdge,VerticesOnGeomEdge,sizeof(VertexOnGeom)*NbVerticesOnGeomEdge);
    4248                 if (NbVertexOnBThEdge)
    4249                  memcpy(newVertexOnBThEdge,VertexOnBThEdge,sizeof(VertexOnEdge)*NbVertexOnBThEdge);
    4250                 Edge *newedges = new Edge [newnbe];
    4251                 //  memcpy(newedges,edges,sizeof(Edge)*nbe);
    4252                 SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv);
    4253                 long k=nbv;
    4254                 long kk=0;
    4255                 long kvb = NbVertexOnBThEdge;
    4256                 long kvg = NbVerticesOnGeomEdge;
    4257                 long ie =0;
    4258                 Edge ** edgesGtoB=0;
    4259                 if (withBackground)
    4260                  edgesGtoB= BTh.MakeGeomEdgeToEdge();
    4261                 long ferr=0;
    4262                 for (i=0;i<nbe;i++)
    4263                  newedges[ie].GeomEdgeHook=0;
    4264 
    4265                 for (i=0;i<nbe;i++)
    4266                   {
    4267                         GeomEdge *ong =  edges[i].GeomEdgeHook;
    4268 
    4269                         newedges[ie]=edges[i];
    4270                         newedges[ie].adj[0]=newedges+(edges[i].adj[0]-edges) ;
    4271                         newedges[ie].adj[1]=newedges + ie +1;
    4272                         R2 A = edges[i][0],B = edges[i][1];
    4273 
    4274                         kk += (i == edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1])));
    4275                         if (ong) // a geometrical edges
    4276                           {
    4277                                 if (withBackground){
    4278                                         // walk on back ground mesh
    4279                                         //  newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge);
    4280                                         // a faire -- difficile
    4281                                         // the first PB is to now a background edge between the 2 vertices
    4282                                         if (!edgesGtoB){
    4283                                                 _error_("!edgesGtoB");
    4284                                         }
    4285                                         ong= ProjectOnCurve(*edgesGtoB[Gh.GetId(edges[i].GeomEdgeHook)],
    4286                                                                 edges[i][0],edges[i][1],0.5,vertices[k],
    4287                                                                 newVertexOnBThEdge[kvb],
    4288                                                                 newVerticesOnGeomEdge[kvg++]);
    4289                                         vertices[k].ReferenceNumber= edges[i].ReferenceNumber;
    4290                                         vertices[k].DirOfSearch =   NoDirOfSearch;       
    4291                                         ;
    4292                                         // get the Info on background mesh
    4293                                         double s =        newVertexOnBThEdge[kvb];
    4294                                         BamgVertex &  bv0  = newVertexOnBThEdge[kvb][0];
    4295                                         BamgVertex &  bv1  = newVertexOnBThEdge[kvb][1];
    4296                                         // compute the metrix of the new points
    4297                                         vertices[k].m =  Metric(1-s,bv0,s,bv1);
    4298                                         kvb++;
    4299                                   }
    4300                                 else
    4301                                   {
    4302                                         ong=Gh.ProjectOnCurve(edges[i],
    4303                                                                 0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);
    4304                                         // vertices[k].i = R2ToI2( vertices[k].r);
    4305                                         vertices[k].ReferenceNumber = edges[i].ReferenceNumber;
    4306                                         vertices[k].DirOfSearch = NoDirOfSearch;
    4307                                         vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);             
    4308                                   } 
    4309                           }
    4310                         else // straigth line edge ---
    4311                           {
    4312                                 vertices[k].r = ((R2) edges[i][0] + (R2)  edges[i][1] )*0.5;
    4313                                 vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);
    4314                                 vertices[k].GeomEdgeHook = 0;
    4315                           }
    4316                         //vertices[k].i = R2ToI2( vertices[k].r);
    4317                         R2 AB =  vertices[k].r;
    4318                         R2 AA = (A+AB)*0.5;
    4319                         R2 BB = (AB+B)*0.5;
    4320                         vertices[k].ReferenceNumber = edges[i].ReferenceNumber;
    4321                         vertices[k].DirOfSearch = NoDirOfSearch;
    4322 
    4323                         newedges[ie].GeomEdgeHook = Gh.Containing(AA,ong);
    4324                         newedges[ie++].v[1]=vertices+k;
    4325 
    4326                         newedges[ie]=edges[i];
    4327                         newedges[ie].adj[0]=newedges + ie -1;
    4328                         newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;
    4329                         newedges[ie].GeomEdgeHook =  Gh.Containing(BB,ong);
    4330                         newedges[ie++].v[0]=vertices+k;
    4331                         k++;
    4332                   }
    4333                 if (edgesGtoB) delete [] edgesGtoB;
    4334                 edgesGtoB=0;
    4335 
    4336                 newnbv=k;
    4337                 newNbVerticesOnGeomEdge=kvg;
    4338                 if (newnbv> maxnbv) goto Error;// bug
    4339 
    4340                 nbv = k;
    4341 
    4342                 kedge = new long[3*nbt+1];
    4343                 ksplitarray = new long[nbt+1];
    4344                 ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0]
    4345 
    4346                 for (i=0;i<3*nbt;i++)
    4347                  kedge[i]=-1;
    4348 
    4349                 // 
    4350 
    4351                 for (i=0;i<nbt;i++) {
    4352                         Triangle & t = triangles[i];
    4353                         if (!t.link){
    4354                                 _error_("!t.link");
    4355                         }
    4356                         for(int j=0;j<3;j++)
    4357                           {
    4358                                 const AdjacentTriangle ta = t.Adj(j);
    4359                                 const Triangle & tt = ta;
    4360                                 if (&tt >= lastT)
    4361                                  t.SetAdj2(j,0,0);// unset adj
    4362                                 const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
    4363                                 const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
    4364                                 long  ke =edge4->SortAndFind(GetId(v0),GetId(v1));
    4365                                 if (ke>0)
    4366                                   {
    4367                                         long ii = GetId(tt);
    4368                                         int  jj = ta;
    4369                                         long ks = ke + nbvold;
    4370                                         kedge[3*i+j] = ks;
    4371                                         if (ii<nbt) // good triangle
    4372                                          kedge[3*ii+jj] = ks;
    4373                                         BamgVertex &A=vertices[ks];
    4374                                         double aa,bb,cc,dd;
    4375                                         if ((dd=Area2(v0.r,v1.r,A.r)) >=0){
    4376                                                 // warning PB roundoff error
    4377                                                 if (t.link && ( (aa=Area2( A.r    , t[1].r , t[2].r )) < 0.0
    4378                                                                                 ||   (bb=Area2( t[0].r , A.r    , t[2].r )) < 0.0 
    4379                                                                                 ||   (cc=Area2( t[0].r , t[1].r , A.r    )) < 0.0)){
    4380                                                         _printf_(ke + nbvold << " not in triangle " << i << " In= " << !!t.link << " " << aa << " " << bb << " " << cc << " " << dd << "\n");
    4381                                                         _error_("Number of triangles with P2 interpolation Problem");
    4382                                                 }
    4383                                         }
    4384                                         else {
    4385                                                 if (tt.link && ( (aa=Area2( A.r     , tt[1].r , tt[2].r )) < 0
    4386                                                                                 ||   (bb=Area2( tt[0].r , A.r     , tt[2].r )) < 0
    4387                                                                                 ||   (cc=Area2( tt[0].r , tt[1].r , A.r     )) < 0)){
    4388                                                         _printf_(ke + nbvold << " not in triangle " << ii << " In= " << !!tt.link << " " << aa << " " << bb << " " << cc << " " << dd << "\n");
    4389                                                         _error_("Number of triangles with P2 interpolation Problem");
    4390                                                 }
    4391                                         }
    4392                                   }
    4393                           }
    4394                 }
    4395 
    4396                 for (i=0;i<nbt;i++){
    4397                         ksplit[i]=1; // no split by default
    4398                         const Triangle & t = triangles[ i];
    4399                         int nbsplitedge =0;
    4400                         int nbinvisible =0;
    4401                         int invisibleedge=0;
    4402                         int kkk[3];     
    4403                         for (int j=0;j<3;j++)
    4404                           {
    4405                                 if (t.Hidden(j)) invisibleedge=j,nbinvisible++;
    4406 
    4407                                 const AdjacentTriangle ta = t.Adj(j);
    4408                                 const Triangle & tt = ta;
    4409 
    4410                                 const BamgVertex & v0 = t[VerticesOfTriangularEdge[j][0]];
    4411                                 const BamgVertex & v1 = t[VerticesOfTriangularEdge[j][1]];
    4412                                 if ( kedge[3*i+j] < 0)
    4413                                   {
    4414                                         long  ke =edge4->SortAndFind(GetId(v0),GetId(v1));
    4415                                         if (ke<0) // new
    4416                                           {
    4417                                                 if (&tt) // internal triangles all the boundary
    4418                                                   { // new internal edges
    4419                                                         long ii = GetId(tt);
    4420                                                         int  jj = ta;
    4421 
    4422                                                         kedge[3*i+j]=k;// save the vertex number
    4423                                                         kedge[3*ii+jj]=k;
    4424                                                         if (k<maxnbv)
    4425                                                           {
    4426                                                                 vertices[k].r = ((R2) v0+(R2) v1 )/2;
    4427                                                                 //vertices[k].i = R2ToI2( vertices[k].r);
    4428                                                                 vertices[k].ReferenceNumber=0;
    4429                                                                 vertices[k].DirOfSearch =NoDirOfSearch;
    4430                                                                 vertices[k].m =  Metric(0.5,v0,0.5,v1);
    4431                                                           }
    4432                                                         k++;
    4433                                                         kkk[nbsplitedge++]=j;                 
    4434                                                   } // tt
    4435                                                 else
    4436                                                  _error_("Bug...");
    4437                                           } // ke<0           
    4438                                         else
    4439                                           { // ke >=0
    4440                                                 kedge[3*i+j]=nbvold+ke;
    4441                                                 kkk[nbsplitedge++]=j;// previously splited
    4442                                           }
    4443                                   }
    4444                                 else
    4445                                  kkk[nbsplitedge++]=j;// previously splited
    4446 
    4447                           }
    4448                         if (nbinvisible>=2){
    4449                                 _error_("nbinvisible>=2");
    4450                         }
    4451                         switch (nbsplitedge) {
    4452                                 case 0: ksplit[i]=10; newnbt++; break;   // nosplit
    4453                                 case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2
    4454                                 case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3
    4455                                 case 3:
    4456                                                   if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4;
    4457                                                   else   ksplit[i]=10*nfortria,newnbt+=nfortria;
    4458                                                   break;
    4459                         }
    4460                         if (ksplit[i]<40){
    4461                                 _error_("ksplit[i]<40");
    4462                         }
    4463                   }
    4464                 //  now do the element split
    4465                 newnbq = 4*nbq;
    4466                 nbv = k;
    4467                 kkk = nbt;
    4468                 ksplit[-1] = nbt;
    4469                 // look on  old true  triangles
    4470 
    4471                 for (i=0;i<nbtsave;i++){
    4472                         int  nbmkadj=0;
    4473                         long mkadj [100];
    4474                         mkadj[0]=i;
    4475                         long kk=ksplit[i]/10;
    4476                         int  ke=(int) (ksplit[i]%10);
    4477                         if (kk>=7 || kk<=0){
    4478                                 _error_("kk>=7 || kk<=0");
    4479                         }
    4480 
    4481                         // def the numbering   k (edge) i vertex
    4482                         int k0 = ke;
    4483                         int k1 = NextEdge[k0];
    4484                         int k2 = PreviousEdge[k0];
    4485                         int i0 = OppositeVertex[k0];
    4486                         int i1 = OppositeVertex[k1];
    4487                         int i2 = OppositeVertex[k2];
    4488 
    4489                         Triangle &t0=triangles[i];
    4490                         BamgVertex * v0=t0(i0);           
    4491                         BamgVertex * v1=t0(i1);           
    4492                         BamgVertex * v2=t0(i2);
    4493 
    4494                         if (nbmkadj>=10){
    4495                                 _error_("nbmkadj>=10");
    4496                         }
    4497                         // --------------------------
    4498                         AdjacentTriangle ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2));
    4499                         // save the flag Hidden
    4500                         int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)};
    4501                         // un set all adj -- save Hidden flag --
    4502                         t0.SetAdj2(0,0,hid[0]);
    4503                         t0.SetAdj2(1,0,hid[1]);
    4504                         t0.SetAdj2(2,0,hid[2]);
    4505                         // --  remake
    4506                         switch  (kk) {
    4507                                 case 1: break;// nothing
    4508                                 case 2: //
    4509                                                   {
    4510                                                         Triangle &t1=triangles[kkk++];
    4511                                                         t1=t0;
    4512                                                         if (kedge[3*i+i0]<0){
    4513                                                                 _error_("kedge[3*i+i0]<0");
    4514                                                         }
    4515                                                         BamgVertex * v3 = vertices + kedge[3*i+k0];
    4516 
    4517                                                         t0(i2) = v3;
    4518                                                         t1(i1) = v3;
    4519                                                         t0.SetAllFlag(k2,0);
    4520                                                         t1.SetAllFlag(k1,0);
    4521                                                   }
    4522                                                 break;
    4523                                 case 3: //
    4524                                                   {
    4525                                                         Triangle &t1=triangles[kkk++];
    4526                                                         Triangle &t2=triangles[kkk++];
    4527                                                         t2=t1=t0;
    4528                                                         if (kedge[3*i+k1]<0){
    4529                                                                 _error_("kedge[3*i+k1]<0");
    4530                                                         }
    4531                                                         if (kedge[3*i+k2]<0){
    4532                                                                 _error_("kedge[3*i+k2]<0");
    4533                                                         }
    4534 
    4535                                                         BamgVertex * v01 = vertices + kedge[3*i+k2];
    4536                                                         BamgVertex * v02 = vertices + kedge[3*i+k1];
    4537                                                         t0(i1) = v01;
    4538                                                         t0(i2) = v02;
    4539                                                         t1(i2) = v02;
    4540                                                         t1(i0) = v01;
    4541                                                         t2(i0) = v02;
    4542                                                         t0.SetAllFlag(k0,0);
    4543                                                         t1.SetAllFlag(k1,0);
    4544                                                         t1.SetAllFlag(k0,0);
    4545                                                         t2.SetAllFlag(k2,0);
    4546                                                   }
    4547                                                 break;
    4548                                 case 4: //
    4549                                 case 6: // split in 4
    4550                                                   {
    4551                                                         Triangle &t1=triangles[kkk++];
    4552                                                         Triangle &t2=triangles[kkk++];
    4553                                                         Triangle &t3=triangles[kkk++];
    4554                                                         t3=t2=t1=t0;
    4555                                                         if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){
    4556                                                                 _error_("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0");
    4557                                                         }
    4558                                                         BamgVertex * v12 = vertices + kedge[3*i+k0];
    4559                                                         BamgVertex * v02 = vertices + kedge[3*i+k1];
    4560                                                         BamgVertex * v01 = vertices + kedge[3*i+k2];
    4561                                                         t0(i1) = v01;
    4562                                                         t0(i2) = v02;
    4563                                                         t0.SetAllFlag(k0,hid[k0]);
    4564 
    4565                                                         t1(i0) = v01;
    4566                                                         t1(i2) = v12;
    4567                                                         t0.SetAllFlag(k1,hid[k1]);
    4568 
    4569                                                         t2(i0) = v02;
    4570                                                         t2(i1) = v12;
    4571                                                         t2.SetAllFlag(k2,hid[k2]);
    4572 
    4573                                                         t3(i0) = v12;
    4574                                                         t3(i1) = v02;
    4575                                                         t3(i2) = v01;
    4576 
    4577                                                         t3.SetAllFlag(0,hid[0]);           
    4578                                                         t3.SetAllFlag(1,hid[1]);           
    4579                                                         t3.SetAllFlag(2,hid[2]);
    4580 
    4581                                                         if ( kk == 6)
    4582                                                           {
    4583 
    4584                                                                 Triangle &t4=triangles[kkk++];
    4585                                                                 Triangle &t5=triangles[kkk++];
    4586 
    4587                                                                 t4 = t3;
    4588                                                                 t5 = t3;
    4589 
    4590                                                                 t0.SetHidden(k0);
    4591                                                                 t1.SetHidden(k1);
    4592                                                                 t2.SetHidden(k2);
    4593                                                                 t3.SetHidden(0);
    4594                                                                 t4.SetHidden(1);
    4595                                                                 t5.SetHidden(2);
    4596 
    4597                                                                 if (nbv < maxnbv )
    4598                                                                   {
    4599                                                                         vertices[nbv].r = ((R2) *v01 + (R2) *v12  + (R2) *v02 ) / 3.0;
    4600                                                                         vertices[nbv].ReferenceNumber =0;
    4601                                                                         vertices[nbv].DirOfSearch =NoDirOfSearch;
    4602                                                                         //vertices[nbv].i = R2ToI2(vertices[nbv].r);
    4603                                                                         double a3[]={1./3.,1./3.,1./3.};
    4604                                                                         vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m);
    4605                                                                         BamgVertex * vc =  vertices +nbv++;
    4606                                                                         t3(i0) = vc;
    4607                                                                         t4(i1) = vc;
    4608                                                                         t5(i2) = vc;
    4609 
    4610                                                                   }
    4611                                                                 else
    4612                                                                  goto Error;
    4613                                                           }
    4614 
    4615                                                   }
    4616                                                 break;         
    4617                         }
    4618 
    4619                         // save all the new triangles
    4620                         mkadj[nbmkadj++]=i;
    4621                         long jj;
    4622                         if (t0.link)
    4623                          for (jj=nbt;jj<kkk;jj++)
    4624                                 {
    4625                                  triangles[jj].link=t0.link;
    4626                                  t0.link= triangles+jj;
    4627                                  mkadj[nbmkadj++]=jj;
    4628                                 }
    4629                         if (nbmkadj>13){// 13 = 6 + 4 +
    4630                                 _error_("nbmkadj>13");
    4631                         }
    4632 
    4633                         if (kk==6)  newnbq+=3;
    4634                         for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;
    4635                         ksplit[i]= nbt; // save last adresse of the new triangles
    4636                         kkk = nbt;
    4637                   }
    4638 
    4639                 for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.;
    4640 
    4641                 if(withBackground)
    4642                  for (i=0;i<BTh.nbv;i++)
    4643                   BTh.vertices[i].m =  BTh.vertices[i].m*2.;
    4644 
    4645                 ret = 2;
    4646                 if (nbt>= maxnbt) goto Error; // bug
    4647                 if (nbv>= maxnbv) goto Error; // bug
    4648                 // generation of the new triangles
    4649 
    4650                 SetIntCoor("In SplitElement");
    4651 
    4652                 CreateSingleVertexToTriangleConnectivity();
    4653                 if(withBackground)
    4654                  BTh.CreateSingleVertexToTriangleConnectivity();
    4655 
    4656                 delete [] edges;
    4657                 edges = newedges;
    4658                 nbe = newnbe;
    4659                 nbq = newnbq;
    4660 
    4661                 for (i=0;i<nbsubdomains;i++)
    4662                   {
    4663                         long k = subdomains[i].edge- edges;
    4664                         subdomains[i].edge =  edges+2*k; // spilt all edge in 2
    4665                   }
    4666 
    4667                 if (ksplitarray) delete [] ksplitarray;
    4668                 if (kedge) delete [] kedge;
    4669                 if (edge4) delete edge4;
    4670                 if (VerticesOnGeomEdge) delete [] VerticesOnGeomEdge;
    4671                 VerticesOnGeomEdge= newVerticesOnGeomEdge;
    4672                 if(VertexOnBThEdge) delete []  VertexOnBThEdge;
    4673                 VertexOnBThEdge = newVertexOnBThEdge;
    4674                 NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
    4675                 NbVertexOnBThEdge=newNbVertexOnBThEdge;
    4676                 //  CreateSingleVertexToTriangleConnectivity();
    4677 
    4678                 ReconstructExistingMesh();
    4679 
    4680                 if (verbose>2){
    4681                         _printf_("   number of quadrilaterals    = " << nbq << "\n");
    4682                         _printf_("   number of triangles         = " << nbt-nbtout- nbq*2 << "\n");
    4683                         _printf_("   number of outside triangles = " << nbtout << "\n");
    4684                 }
    4685 
    4686                 return 0; //ok
    4687 
    4688 Error:
    4689                 nbv = nbvold;
    4690                 nbt = nbtold;
    4691                 nbtout = nbtoutold;
    4692                 // cleaning memory ---
    4693                 delete [] newedges;
    4694                 if (ksplitarray) delete [] ksplitarray;
    4695                 if (kedge) delete [] kedge;
    4696                 if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge;
    4697                 if (edge4) delete edge4;
    4698                 if(newVertexOnBThEdge) delete []  newVertexOnBThEdge;
    4699 
    4700                 return ret; // ok
    4701         }
    4702         /*}}}*/
    47033987/*FUNCTION Mesh::SplitInternalEdgeWithBorderVertices{{{*/
    47043988long  Mesh::SplitInternalEdgeWithBorderVertices(){
    … …  
    48884172        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)*/
    48894173
    4890         long int verbose=0;
    48914174        register Triangle *t0,*t;
    48924175        register long k=0, num;   
    … …  
    49374220        //Vertices
    49384221        if(verbose) _printf_("Reading vertices (" << nbv << ")\n");
    4939         for (i=0;i<nbv;i++){
     4222        for(i=0;i<nbv;i++){
    49404223                vertices[i].r.x=x[i];
    49414224                vertices[i].r.y=y[i];
    … …  
    54154698
    54164699                                /*Initialize variables*/
    5417                                 double Lstep=0,Lcurve=0;    // step between two points   (phase==1)
     4700                                double Lstep=0;             // step between two points   (phase==1)
    54184701                                long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
    54194702
    … …  
    55484831                                                long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
    55494832                                                Lstep = L/NbSegOnCurve;
    5550                                                 Lcurve = L;
    55514833                                                NbCreatePointOnCurve = NbSegOnCurve-1;
    55524834                                                NbOfNewEdge += NbSegOnCurve;
    … …  
    56794961                                _error_("!v1 || !v2");
    56804962                        }
    5681                         Icoor2 detss = 0,l=0,ks;
    5682                         while ((ks=SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
     4963                        Icoor2 detss = 0,l=0;
     4964                        while ((SwapForForcingEdge(  va,  vb, tc, detss, det1,det2,NbSwap)))
    56834965                         if(l++ > 10000000) {
    56844966                                 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");
    56854967                         }
    56864968                        BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
    5687                         if (( aa == &a ) && (bb == &b) ||  (bb ==  &a ) && (aa == &b)) {
     4969                        if (((aa == &a ) && (bb == &b)) ||((bb ==  &a ) && (aa == &b))){
    56884970                                tc.SetLock();
    56894971                                a.Optim(1,0);
    … …  
    58275109                        detsb=-detsb;
    58285110
    5829                         if (ToSwap)
     5111                        if(ToSwap){
    58305112                         if (dets2 < 0) {// haut
    58315113                                 dets1 = (ToSwap ? dets1 : detsa) ;
    58325114                                 detsa = dets2;
    58335115                                 tt1 =  Previous(tt2) ;}
    5834                          else if (dets2 > 0){// bas
     5116                         else if(dets2 > 0){// bas
    58355117                                 dets1 = (ToSwap ? dets1 : detsb) ;
    58365118                                 detsb =  dets2;
    … …  
    58405122                                 tt1 = Next(tt2);
    58415123                                 ret =0;}
     5124                        }
    58425125
    58435126                }
Note: See TracChangeset for help on using the changeset viewer.