Changeset 5148


Ignore:
Timestamp:
08/11/10 10:38:41 (15 years ago)
Author:
Mathieu Morlighem
Message:

as usual

Location:
issm/trunk/src/c
Files:
22 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/modules/Bamgx/Bamgx.cpp

    r5146 r5148  
    184184                //display info
    185185                if(verbosity>0) {
    186                         if (Th.nbt-Th.nbtout-Th.NbOfQuad*2){
    187                                 printf("   new number of triangles = %i\n",(Th.nbt-Th.nbtout-Th.NbOfQuad*2));
    188                         }
    189                         if (Th.NbOfQuad ){
    190                                 printf("   new number of quads = %i\n",Th.NbOfQuad);
     186                        if (Th.nbt-Th.nbtout-Th.nbq*2){
     187                                printf("   new number of triangles = %i\n",(Th.nbt-Th.nbtout-Th.nbq*2));
     188                        }
     189                        if (Th.nbq ){
     190                                printf("   new number of quads = %i\n",Th.nbq);
    191191                        }
    192192                }
  • issm/trunk/src/c/objects/Bamg/AdjacentTriangle.cpp

    r5143 r5148  
    1616        /*FUNCTION AdjacentTriangle::Locked {{{1*/
    1717        int  AdjacentTriangle::Locked() const {
    18                 return t->TriaAdjSharedEdge[a] & 4;
     18                return t->AdjEdgeNumber[a] & 4;
    1919        }
    2020        /*}}}*/
    2121        /*FUNCTION AdjacentTriangle::MarkUnSwap {{{1*/
    2222        int  AdjacentTriangle::MarkUnSwap() const {
    23                 return t->TriaAdjSharedEdge[a] & 8;
     23                return t->AdjEdgeNumber[a] & 8;
    2424        }
    2525        /*}}}*/
     
    2727        int  AdjacentTriangle::GetAllFlag_UnSwap() const {
    2828                // take all flag except MarkUnSwap
    29                 return t->TriaAdjSharedEdge[a] & 1012;
     29                return t->AdjEdgeNumber[a] & 1012;
    3030        }
    3131        /*}}}*/
     
    4242        /*FUNCTION AdjacentTriangle::EdgeVertex {{{1*/
    4343        BamgVertex* AdjacentTriangle::EdgeVertex(const int & i) const {
    44                 return t->TriaVertices[VerticesOfTriangularEdge[a][i]];
     44                return t->vertices[VerticesOfTriangularEdge[a][i]];
    4545        }
    4646        /*}}}*/
    4747        /*FUNCTION AdjacentTriangle::OppositeVertex {{{1*/
    4848        BamgVertex* AdjacentTriangle::OppositeVertex() const {
    49                 return t->TriaVertices[bamg::OppositeVertex[a]];
     49                return t->vertices[bamg::OppositeVertex[a]];
    5050        }
    5151        /*}}}*/
     
    6464                //set Adjacent Triangle of a triangle
    6565                if(t) {
    66                         t->TriaAdjTriangles[a]=ta.t;
    67                         t->TriaAdjSharedEdge[a]=ta.a|l;
     66                        t->adj[a]=ta.t;
     67                        t->AdjEdgeNumber[a]=ta.a|l;
    6868                }
    6969                if(ta.t) {
    70                         ta.t->TriaAdjTriangles[ta.a] = t ;
    71                         ta.t->TriaAdjSharedEdge[ta.a] = a| l ;
     70                        ta.t->adj[ta.a] = t ;
     71                        ta.t->AdjEdgeNumber[ta.a] = a| l ;
    7272                }
    7373        }
  • issm/trunk/src/c/objects/Bamg/BamgVertex.h

    r5130 r5148  
    2828
    2929                        union {
    30                                 Triangle     *t;                    // one triangle which is containing the vertex
     30                                Triangle     *t;                      // one triangle which is containing the vertex
    3131                                long          color;
    32                                 BamgVertex   *to;                   // used in geometry BamgVertex to know the Mesh BamgVertex associated
    33                                 VertexOnGeom *onGeometry;           // if vint == 8; // set with Mesh::SetVertexFieldOn()
    34                                 BamgVertex   *onBackgroundVertex;   // if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh()
    35                                 VertexOnEdge *onBackgroundEdge;     // if vint == 32 on Background edge
     32                                BamgVertex   *to;                     // used in geometry BamgVertex to know the Mesh Vertex associated
     33                                VertexOnGeom *GeometricalEdgeHook;    // if vint == 8; // set with Mesh::SetVertexFieldOn()
     34                                BamgVertex   *BackgroundVertexHook;   // if vint == 16 on Background vertex Mesh::SetVertexFieldOnBTh()
     35                                VertexOnEdge *BackgroundEdgeHook;     // if vint == 32 on Background edge
    3636                        };
    3737
     
    4646                        void   MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts);
    4747                        void   Echo();
    48                         int    ref() const { return ReferenceNumber;}
     48                        int    GetReferenceNumber() const { return ReferenceNumber;}
    4949                        long   Optim(int =1,int =0);
    5050
  • issm/trunk/src/c/objects/Bamg/Direction.cpp

    r5130 r5148  
    2323
    2424        /*Methods*/
    25         /*FUNCTION Direction::sens{{{1*/
    26         int Direction::sens(Icoor1 i,Icoor1 j) {
     25        /*FUNCTION Direction::direction{{{1*/
     26        int Direction::direction(Icoor1 i,Icoor1 j) {
    2727                int r =1;
    2828                if (dir!= MaxICoor) {
  • issm/trunk/src/c/objects/Bamg/Direction.h

    r5130 r5148  
    1515                        Direction();
    1616                        Direction(Icoor1 i,Icoor1 j);
    17                         int sens(Icoor1 i,Icoor1 j);
     17                        int direction(Icoor1 i,Icoor1 j);
    1818        };
    1919}
  • issm/trunk/src/c/objects/Bamg/Edge.cpp

    r5130 r5148  
    1818                v[0] = ThNew.vertices + Th.Number(v[0]);   
    1919                v[1] = ThNew.vertices + Th.Number(v[1]);
    20                 if (onGeometry)
    21                  onGeometry =  ThNew.Gh.edges+Th.Gh.Number(onGeometry);
     20                if (GeometricalEdgeHook)
     21                 GeometricalEdgeHook =  ThNew.Gh.edges+Th.Gh.Number(GeometricalEdgeHook);
    2222                if (adj[0]) adj[0] =   ThNew.edges +   Th.Number(adj[0]);
    2323                if (adj[1]) adj[1] =   ThNew.edges +   Th.Number(adj[1]);
     
    2828                printf("Edge:\n");
    2929                printf("   pointers towards two vertices: %p %p\n",v[0],v[1]);
    30                 printf("   ref = %i\n",ref);
    31                 printf("   onGeometry = %p\n",onGeometry);
     30                printf("   ReferenceNumber = %i\n",ReferenceNumber);
     31                printf("   GeometricalEdgeHook = %p\n",GeometricalEdgeHook);
    3232                printf("   two adjacent edges on the same curve: %p %p\n",adj[0],adj[1]);
    3333        }
  • issm/trunk/src/c/objects/Bamg/Edge.h

    r5143 r5148  
    1717                public:
    1818                        BamgVertex      *v[2];
    19                         long             ref;
    20                         GeometricalEdge *onGeometry;
     19                        long             ReferenceNumber;
     20                        GeometricalEdge *GeometricalEdgeHook;
    2121                        Edge            *adj[2];       // the 2 adj edges if on the same curve
    2222
  • issm/trunk/src/c/objects/Bamg/GeometricalEdge.h

    r5146 r5148  
    1414                public:
    1515                        GeometricalVertex *v[2];
    16                         long               ref;
     16                        long               ReferenceNumber;
    1717                        long               CurveNumber;
    1818                        R2                 tg[2];         // the 2 tangentes (tg[0] =0 => no continuity)
     
    2929                        R2     F(double theta) const ; // parametrization of the curve edge
    3030                        double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente
    31                         int    Tg(int i) const{return i==0 ? TgA() : TgB();}
    32                         int    Cracked() const{return type & 1;}
    33                         int    TgA()const{return type & 4; }
    34                         int    TgB()const {return type & 8;}
    35                         int    Mark()const {return type & 16;}
    36                         int    Required() {return type & 64;}
     31                        int    Tg(int i) const  {return i==0 ? TgA() : TgB(); }
     32                        int    Cracked() const  {return type &1;  }
     33                        int    TgA()     const  {return type &4; }
     34                        int    TgB()     const  {return type &8;  }
     35                        int    Mark()    const  {return type &16; }
     36                        int    Required()       {return type &64; }
    3737                        void   SetCracked()     { type |= 1;}
    3838                        void   SetTgA()         { type |=4; }
  • issm/trunk/src/c/objects/Bamg/GeometricalSubDomain.h

    r3913 r5148  
    1313                public:
    1414                        GeometricalEdge *edge;
    15                         int sens; // -1 or 1
    16                         long ref;
     15                        int              direction;  // -1 or 1
     16                        long             ReferenceNumber;
    1717
    1818                        //Methods
  • issm/trunk/src/c/objects/Bamg/Geometry.cpp

    r5146 r5148  
    3535                vertices = nbv ? new GeometricalVertex[nbv] : NULL;
    3636                edges = nbe ? new GeometricalEdge[nbe]:NULL;
    37                 curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;
    38                 subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL;
     37                curves= nbcurves ? new Curve[nbcurves]:NULL;
     38                subdomains = nbsubdomains ? new GeometricalSubDomain[nbsubdomains]:NULL;
    3939                for (i=0;i<nbe;i++)
    4040                 edges[i].Set(Gh.edges[i],Gh,*this);
    41                 for (i=0;i<NbOfCurves;i++)
     41                for (i=0;i<nbcurves;i++)
    4242                 curves[i].Set(Gh.curves[i],Gh,*this);
    43                 for (i=0;i<NbSubDomains;i++)
     43                for (i=0;i<nbsubdomains;i++)
    4444                 subdomains[i].Set(Gh.subdomains[i],Gh,*this);
    4545        }
     
    5151                if(vertices)  delete [] vertices;   vertices=0;
    5252                if(edges)     delete [] edges;      edges=0;
    53                 if(triangles) delete [] triangles;  triangles=0;
    5453                if(quadtree)  delete  quadtree;     quadtree=0;
    55                 if(curves)    delete  []curves;     curves=0;NbOfCurves=0;
     54                if(curves)    delete  []curves;     curves=0;nbcurves=0;
    5655                if(subdomains) delete [] subdomains;subdomains=0;
    5756                Init();
     
    6665                nbv=0;
    6766                nbe=0;
    68                 NbOfCurves=0;
     67                nbcurves=0;
    6968
    7069                double Hmin = HUGE_VAL;// the infinie value
     
    140139                                edges[i].v[0]= vertices + i1;     //pointer toward vertex i1 (=&vertices[i1])
    141140                                edges[i].v[1]= vertices + i2;     //pointer toward vertex i2
    142                                 edges[i].ref=(long)bamggeom->Edges[i*3+2];
     141                                edges[i].ReferenceNumber=(long)bamggeom->Edges[i*3+2];
    143142
    144143                                //get length of edge
     
    155154
    156155                                //Cracked?
    157                                 if (edges[i].ref!=1) edges[i].SetCracked();
     156                                if (edges[i].ReferenceNumber!=1) edges[i].SetCracked();
    158157
    159158                                //prepare metric
     
    263262                        if(verbose>5) printf("      processing SubDomains\n");
    264263                        if (bamggeom->SubDomainsSize[1]!=4) ISSMERROR("SubDomains should have 4 columns");
    265                         NbSubDomains=bamggeom->SubDomainsSize[0];
    266                         subdomains = new GeometricalSubDomain[NbSubDomains];
    267                         for (i=0;i<NbSubDomains;i++){
     264                        nbsubdomains=bamggeom->SubDomainsSize[0];
     265                        subdomains = new GeometricalSubDomain[nbsubdomains];
     266                        for (i=0;i<nbsubdomains;i++){
    268267                                i0=(int)bamggeom->SubDomains[i*4+0];
    269268                                i1=(int)bamggeom->SubDomains[i*4+1];
     
    273272                                if (i1>nbe || i1<=0) ISSMERROR("Bad Subdomain definition: second number should in [1 %i] (edge number)",nbe);
    274273                                subdomains[i].edge=edges + (i1-1);
    275                                 subdomains[i].sens = (int) i2;
    276                                 subdomains[i].ref = i3;
     274                                subdomains[i].direction = (int) i2;
     275                                subdomains[i].ReferenceNumber = i3;
    277276                        }
    278277                }
     
    304303                                bamggeom->Vertices[i*3+0]=vertices[i].r.x;
    305304                                bamggeom->Vertices[i*3+1]=vertices[i].r.y;
    306                                 bamggeom->Vertices[i*3+2]=vertices[i].ref();
     305                                bamggeom->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
    307306
    308307                                //update counters
     
    320319                                bamggeom->Edges[i*3+0]=Number(edges[i][0])+1; //back to Matlab indexing
    321320                                bamggeom->Edges[i*3+1]=Number(edges[i][1])+1; //back to Matlab indexing
    322                                 bamggeom->Edges[i*3+2]=(double)edges[i].ref;
     321                                bamggeom->Edges[i*3+2]=(double)edges[i].ReferenceNumber;
    323322
    324323                                //update counters
     
    363362                /*SubDomains*/
    364363                if(verbose>5) printf("      writing SubDomains\n");
    365                 bamggeom->SubDomainsSize[0]=NbSubDomains;
     364                bamggeom->SubDomainsSize[0]=nbsubdomains;
    366365                bamggeom->SubDomainsSize[1]=4;
    367                 if (NbSubDomains){
    368                         bamggeom->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));
    369                         for (i=0;i<NbSubDomains;i++){
     366                if (nbsubdomains){
     367                        bamggeom->SubDomains=(double*)xmalloc(4*nbsubdomains*sizeof(double));
     368                        for (i=0;i<nbsubdomains;i++){
    370369                                bamggeom->SubDomains[4*i+0]=2;
    371370                                bamggeom->SubDomains[4*i+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing
    372                                 bamggeom->SubDomains[4*i+2]=subdomains[i].sens;
    373                                 bamggeom->SubDomains[4*i+3]=subdomains[i].ref;
     371                                bamggeom->SubDomains[4*i+2]=subdomains[i].direction;
     372                                bamggeom->SubDomains[4*i+3]=subdomains[i].ReferenceNumber;
    374373                        }
    375374                }
     
    585584                                        vertices[i].SetCorner() ;
    586585                                }
    587                                 // if the ref a changing then is     SetRequired();
     586                                // if the ReferenceNumber a changing then is     SetRequired();
    588587                                if (edges[i1].type != edges[i2].type || edges[i1].Required()){
    589588                                        vertices[i].SetRequired();
    590589                                }
    591                                 if (edges[i1].ref != edges[i2].ref) {
     590                                if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
    592591                                        vertices[i].SetRequired();
    593592                                }
     
    661660                for (int step=0;step<2;step++){
    662661                        for (i=0;i<nbe;i++) edges[i].SetUnMark();
    663                         NbOfCurves = 0;
     662                        nbcurves = 0;
    664663                        long  nbgem=0;
    665664                        for (int level=0;level < 2 && nbgem != nbe;level++)
     
    673672                                                 GeometricalVertex *a=(*e)(k0); // begin
    674673                                                 if(curves){
    675                                                          curves[NbOfCurves].be=e;
    676                                                          curves[NbOfCurves].kb=k0;
     674                                                         curves[nbcurves].be=e;
     675                                                         curves[nbcurves].kb=k0;
    677676                                                 }
    678677                                                 int nee=0;
     
    682681                                                         e->SetMark();
    683682                                                         nbgem++;
    684                                                          e->CurveNumber=NbOfCurves;
     683                                                         e->CurveNumber=nbcurves;
    685684                                                         if(curves) {
    686                                                                  curves[NbOfCurves].ee=e;
    687                                                                  curves[NbOfCurves].ke=k1;
     685                                                                 curves[nbcurves].ee=e;
     686                                                                 curves[nbcurves].ke=k1;
    688687                                                         }
    689688
     
    693692                                                         e = e->Adj[k1]; // next edge
    694693                                                 }
    695                                                  NbOfCurves++;
     694                                                 nbcurves++;
    696695                                                 if(level) a->SetRequired();
    697696                                         }
     
    700699                        ISSMASSERT(nbgem && nbe);
    701700                        if(step==0){
    702                                 curves = new Curve[NbOfCurves];
     701                                curves = new Curve[nbcurves];
    703702                        }
    704703                }
    705                 for(int i=0;i<NbOfCurves ;i++){
     704                for(int i=0;i<nbcurves ;i++){
    706705                        GeometricalEdge * be=curves[i].be, *eqbe=be;
    707706                        //GeometricalEdge * ee=curves[i].ee, *eqee=be;
     
    749748                printf("   nbv  (number of vertices) : %i\n",nbv);
    750749                printf("   nbe  (number of edges)    : %i\n",nbe);
    751                 printf("   NbSubDomains: %i\n",NbSubDomains);
    752                 printf("   NbOfCurves: %i\n",NbOfCurves);
     750                printf("   nbsubdomains: %i\n",nbsubdomains);
     751                printf("   nbcurves: %i\n",nbcurves);
    753752                printf("   vertices: %p\n",vertices);
    754                 printf("   triangles: %p\n",triangles);
    755753                printf("   edges: %p\n",edges);
    756754                printf("   quadtree: %p\n",quadtree);
     
    774772                quadtree=0;
    775773                curves=0;
    776                 triangles=0;
    777774                edges=0;
    778775                vertices=0;
    779                 NbSubDomains=0;
    780                 NbOfCurves=0;
     776                nbsubdomains=0;
     777                nbcurves=0;
    781778                subdomains=0;
    782779                MaxCornerAngle = 10*Pi/180; //default is 10 degres
     
    822819
    823820                s=save_s;
    824                 GeometricalEdge* on=e.onGeometry;
     821                GeometricalEdge* on=e.GeometricalEdgeHook;
    825822                if (!on){
    826823                        ISSMERROR("ProjectOnCurve error message: edge provided should be on geometry");
    827824                }
    828                 if (!e[0].onGeometry ||  !e[1].onGeometry){
     825                if (!e[0].GeometricalEdgeHook ||  !e[1].GeometricalEdgeHook){
    829826                        ISSMERROR("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry");
    830827                }
     
    838835
    839836                //Get geometrical vertices corresponding to v0 and v1
    840                 VertexOnGeom  vg0=*v0.onGeometry,  vg1=*v1.onGeometry;
     837                VertexOnGeom  vg0=*v0.GeometricalEdgeHook,  vg1=*v1.GeometricalEdgeHook;
    841838
    842839                //build two pointers towrad current geometrical edge
     
    846843                R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
    847844                int OppositeSens = (V01,AB)<0;
    848                 int sens0=0,sens1=1;
     845                int direction0=0,direction1=1;
    849846                if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
    850847
     
    854851                const int mxe=100;
    855852                GeometricalEdge* ge[mxe+1];
    856                 int     sensge[mxe+1];
     853                int     directionge[mxe+1];
    857854                double  lge[mxe+1];
    858855                int bge=mxe/2,tge=bge;
    859                 ge[bge] = e.onGeometry;
    860                 sensge[bge]=1;
    861 
    862                 while (eg0!=(GeometricalEdge*)vg0 && (*eg0)(sens0)!=(GeometricalVertex*)vg0){
     856                ge[bge] = e.GeometricalEdgeHook;
     857                directionge[bge]=1;
     858
     859                while (eg0!=(GeometricalEdge*)vg0 && (*eg0)(direction0)!=(GeometricalVertex*)vg0){
    863860                        if (bge<=0) {
    864861                                if(NbTry) {
    865                                         printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n");
     862                                        printf("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve\n");
    866863                                        printf("That bug might come from:\n");
    867864                                        printf(" 1)  a mesh edge  containing more than %i geometrical edges\n",mxe/2);
     
    874871                        }
    875872                        GeometricalEdge* tmpge = eg0;
    876                         ge[--bge] =eg0 = eg0->Adj[sens0];
     873                        ge[--bge] =eg0 = eg0->Adj[direction0];
    877874                        ISSMASSERT(bge>=0 && bge<=mxe);
    878                         sens0 = 1-( sensge[bge] = tmpge->DirAdj[sens0]);
    879                 }
    880                 while (eg1 != (GeometricalEdge*) vg1  &&  (*eg1)(sens1) != (GeometricalVertex*) vg1) {
     875                        direction0 = 1-( directionge[bge] = tmpge->DirAdj[direction0]);
     876                }
     877                while (eg1 != (GeometricalEdge*) vg1  &&  (*eg1)(direction1) != (GeometricalVertex*) vg1) {
    881878                        if(tge>=mxe ) {
    882                                 printf("WARNING: on the class triangles before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n");
     879                                printf("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n");
    883880                                NbTry++;
    884881                                if (NbTry<2) goto retry;
    885                                 printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n");
     882                                printf("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve\n");
    886883                                printf("That bug might come from:\n");
    887884                                printf(" 1)  a mesh edge  contening more than %i geometrical edges\n",mxe/2);
     
    891888                        }
    892889                        GeometricalEdge* tmpge = eg1;
    893                         ge[++tge] =eg1 = eg1->Adj[sens1];
    894                         sensge[tge]= sens1 = 1-tmpge->DirAdj[sens1];
     890                        ge[++tge] =eg1 = eg1->Adj[direction1];
     891                        directionge[tge]= direction1 = 1-tmpge->DirAdj[direction1];
    895892                        ISSMASSERT(tge>=0 && tge<=mxe);
    896893                }
    897894
    898895
    899                 if ((*eg0)(sens0)==(GeometricalVertex*)vg0)
    900                  vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,sens0); //vg0 = absisce
    901 
    902                 if ((*eg1)(sens1)==(GeometricalVertex*)vg1)
    903                  vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,sens1);
     896                if ((*eg0)(direction0)==(GeometricalVertex*)vg0)
     897                 vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce
     898
     899                if ((*eg1)(direction1)==(GeometricalVertex*)vg1)
     900                 vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,direction1);
    904901
    905902                double sg;
     
    916913                        for(i=bge;i<tge;i++){
    917914                                ISSMASSERT(i>=0 && i<=mxe);
    918                                 BB =  (*ge[i])[sensge[i]];
     915                                BB =  (*ge[i])[directionge[i]];
    919916                                lge[i]=ll += Norme2(AA-BB);
    920917                                AA=BB ;}
     
    925922                                on =0;
    926923                                s0 = vg0;
    927                                 s1= sensge[bge];
     924                                s1= directionge[bge];
    928925                                double l0=0,l1;
    929926                                i=bge;
    930927                                while (  (l1=lge[i]) < ls ) {
    931928                                        ISSMASSERT(i>=0 && i<=mxe);
    932                                         i++,s0=1-(s1=sensge[i]),l0=l1;
     929                                        i++,s0=1-(s1=directionge[i]),l0=l1;
    933930                                }
    934931                                on=ge[i];
  • issm/trunk/src/c/objects/Bamg/Geometry.h

    r5130 r5148  
    2323                        long                  nbv;                           // number of vertices
    2424                        long                  nbe;                           // number of edges
    25                         long                  NbSubDomains;
    26                         long                  NbOfCurves;
     25                        long                  nbsubdomains;
     26                        long                  nbcurves;
    2727                        GeometricalVertex    *vertices;
    28                         Triangle             *triangles;
    2928                        GeometricalEdge      *edges;
    3029                        QuadTree             *quadtree;
  • issm/trunk/src/c/objects/Bamg/Mesh.cpp

    r5146 r5148  
    106106                          {
    107107                                vertices[nbv] = Tho.vertices[i];
    108                                 if (!vertices[nbv].ref())
     108                                if (!vertices[nbv].GetReferenceNumber())
    109109                                 vertices[nbv].ReferenceNumber = refv[i];
    110110                                nbv++;
     
    127127                                }
    128128                                triangles[nbt] = Triangle(this,kk[i0],kk[i1],kk[i2]);
    129                                 triangles[nbt].color = Tho.subdomains[reft[i]].ref;
     129                                triangles[nbt].color = Tho.subdomains[reft[i]].ReferenceNumber;
    130130                                nbt++;           
    131131                          }
     
    145145                  ReconstructExistingMesh();
    146146
    147                   if (!NbSubDomains){
    148                           ISSMERROR("NbSubDomains==0");
     147                  if (!nbsubdomains){
     148                          ISSMERROR("nbsubdomains==0");
    149149                  }
    150150                  if (!subdomains[0].head || !subdomains[0].head->link){
     
    168168                  nbt = Th.nbt;
    169169                  nbe = Th.nbe;
    170                   NbSubDomains = Th.NbSubDomains;
     170                  nbsubdomains = Th.nbsubdomains;
    171171                  nbtout = Th.nbtout;
    172                   NbOfQuad =  Th.NbOfQuad ;
    173                   NbOfSwapTriangle =0;
     172                  nbq =  Th.nbq ;
    174173                  NbVerticesOnGeomVertex = Th.NbVerticesOnGeomVertex;
    175174                  if(NbVerticesOnGeomVertex)
     
    198197                  if(nbe)
    199198                        edges = new Edge[nbe];
    200                   if(NbSubDomains)
    201                         subdomains = new SubDomain[NbSubDomains];
     199                  if(nbsubdomains)
     200                        subdomains = new SubDomain[nbsubdomains];
    202201                  pmin = Th.pmin;
    203202                  pmax = Th.pmax;
     
    209208                  for(i=0;i<nbv;i++)
    210209                        vertices[i].Set(Th.vertices[i],Th,*this);
    211                   for(i=0;i<NbSubDomains;i++) 
     210                  for(i=0;i<nbsubdomains;i++) 
    212211                        subdomains[i].Set(Th,i,*this);
    213212                  for (i=0;i<NbVerticesOnGeomVertex;i++)
     
    424423                                i1=(int)bamgmesh->Edges[i*3+0]-1; //-1 for C indexing
    425424                                i2=(int)bamgmesh->Edges[i*3+1]-1; //-1 for C indexing
    426                                 edges[i].ref=(long)bamgmesh->Edges[i*3+2];
     425                                edges[i].ReferenceNumber=(long)bamgmesh->Edges[i*3+2];
    427426                                edges[i].v[0]= vertices +i1;
    428427                                edges[i].v[1]= vertices +i2;
     
    484483                                        ISSMERROR("ReadMesh error: EdgesOnGeometricEdge edge provided (line %i: [%i %i]) is incorrect (must be positive, [0<i<nbe=%i 0<j<Gh.nbe=%i]",i1+1,i+1,j+1,nbe,Gh.nbe);
    485484                                }
    486                                 edges[i].onGeometry=Gh.edges+j;
     485                                edges[i].GeometricalEdgeHook=Gh.edges+j;
    487486                        }
    488487                }
     
    500499                //SubDomain
    501500                if(bamgmesh->SubDomains){
    502                         long i3,head,sens;
     501                        long i3,head,direction;
    503502                        if(verbose>5) printf("      processing SubDomains\n");
    504                         NbSubDomains=bamgmesh->SubDomainsSize[0];
    505                         subdomains = new SubDomain [ NbSubDomains ];
    506                         for (i=0;i<NbSubDomains;i++) {
     503                        nbsubdomains=bamgmesh->SubDomainsSize[0];
     504                        subdomains = new SubDomain [ nbsubdomains ];
     505                        for (i=0;i<nbsubdomains;i++) {
    507506                                i3  =(int)bamgmesh->SubDomains[i*3+0];
    508507                                head=(int)bamgmesh->SubDomains[i*3+1]-1;//C indexing
    509                                 sens=(int)bamgmesh->SubDomains[i*3+2];
     508                                direction=(int)bamgmesh->SubDomains[i*3+2];
    510509                                if (i3!=23) ISSMERROR("Bad Subdomain definition: first number should be 3");
    511510                                if (head<0 || head>=nbt) ISSMERROR("Bad Subdomain definition: head should in [1 %i] (triangle number)",nbt);
     
    581580                                bamgmesh->Vertices[i*3+0]=vertices[i].r.x;
    582581                                bamgmesh->Vertices[i*3+1]=vertices[i].r.y;
    583                                 bamgmesh->Vertices[i*3+2]=vertices[i].ref();
     582                                bamgmesh->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
    584583                        }
    585584                }
     
    595594                                bamgmesh->Edges[i*3+0]=Number(edges[i][0])+1; //back to M indexing
    596595                                bamgmesh->Edges[i*3+1]=Number(edges[i][1])+1; //back to M indexing
    597                                 bamgmesh->Edges[i*3+2]=edges[i].ref;
    598                                 if(edges[i].onGeometry){
     596                                bamgmesh->Edges[i*3+2]=edges[i].ReferenceNumber;
     597                                if(edges[i].GeometricalEdgeHook){
    599598                                        NumIssmSegments++;
    600599                                }
     
    665664                num=0;
    666665                for (i=0;i<nbe;i++){
    667                         if(edges[i].onGeometry){
     666                        if(edges[i].GeometricalEdgeHook){
    668667                                //build segment
    669668                                int i1=Number(edges[i][0]);
     
    677676                                                                bamgmesh->IssmSegments[num*4+1]=Number(edges[i][1])+1; //back to M indexing
    678677                                                                bamgmesh->IssmSegments[num*4+2]=(int)j/3+1;            //back to M indexing
    679                                                                 bamgmesh->IssmSegments[num*4+3]=edges[i].ref;
     678                                                                bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber;
    680679                                                                num+=1;
    681680                                                                stop=true;
     
    686685                                                                bamgmesh->IssmSegments[num*4+1]=Number(edges[i][0])+1; //back to M indexing
    687686                                                                bamgmesh->IssmSegments[num*4+2]=(int)j/3+1;            //back to M indexing
    688                                                                 bamgmesh->IssmSegments[num*4+3]=edges[i].ref;
     687                                                                bamgmesh->IssmSegments[num*4+3]=edges[i].ReferenceNumber;
    689688                                                                num+=1;
    690689                                                                stop=true;
     
    703702                /*Triangles*/
    704703                if(verbose>5) printf("      writing Triangles\n");
    705                 k=nbInT-NbOfQuad*2;
     704                k=nbInT-nbq*2;
    706705                num=0;
    707706                bamgmesh->TrianglesSize[0]=k;
     
    716715                                        bamgmesh->Triangles[num*4+1]=Number(t[1])+1; //back to M indexing
    717716                                        bamgmesh->Triangles[num*4+2]=Number(t[2])+1; //back to M indexing
    718                                         bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ref;
     717                                        bamgmesh->Triangles[num*4+3]=subdomains[reft[i]].ReferenceNumber;
    719718                                        num=num+1;
    720719                                }
     
    724723                /*Quadrilaterals*/
    725724                if(verbose>5) printf("      writing Quadrilaterals\n");
    726                 bamgmesh->QuadrilateralsSize[0]=NbOfQuad;
     725                bamgmesh->QuadrilateralsSize[0]=nbq;
    727726                bamgmesh->QuadrilateralsSize[1]=5;
    728                 if (NbOfQuad){
    729                         bamgmesh->Quadrilaterals=(double*)xmalloc(5*NbOfQuad*sizeof(double));
     727                if (nbq){
     728                        bamgmesh->Quadrilaterals=(double*)xmalloc(5*nbq*sizeof(double));
    730729                        for (i=0;i<nbt;i++){
    731730                                Triangle &t =triangles[i];
     
    738737                                        bamgmesh->Quadrilaterals[i*5+2]=Number(v2)+1; //back to M indexing
    739738                                        bamgmesh->Quadrilaterals[i*5+3]=Number(v3)+1; //back to M indexing
    740                                         bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ref;
     739                                        bamgmesh->Quadrilaterals[i*5+4]=subdomains[reft[i]].ReferenceNumber;
    741740                                }
    742741                        }
     
    745744                /*SubDomains*/
    746745                if(verbose>5) printf("      writing SubDomains\n");
    747                 bamgmesh->SubDomainsSize[0]=NbSubDomains;
     746                bamgmesh->SubDomainsSize[0]=nbsubdomains;
    748747                bamgmesh->SubDomainsSize[1]=4;
    749                 if (NbSubDomains){
    750                         bamgmesh->SubDomains=(double*)xmalloc(4*NbSubDomains*sizeof(double));
    751                         for (i=0;i<NbSubDomains;i++){
     748                if (nbsubdomains){
     749                        bamgmesh->SubDomains=(double*)xmalloc(4*nbsubdomains*sizeof(double));
     750                        for (i=0;i<nbsubdomains;i++){
    752751                                bamgmesh->SubDomains[i*4+0]=3;
    753752                                bamgmesh->SubDomains[i*4+1]=reft[Number(subdomains[i].head)];
    754753                                bamgmesh->SubDomains[i*4+2]=1;
    755                                 bamgmesh->SubDomains[i*4+3]=subdomains[i].ref;
     754                                bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber;
    756755                        }
    757756                }
     
    759758                /*SubDomainsFromGeom*/
    760759                if(verbose>5) printf("      writing SubDomainsFromGeom\n");
    761                 bamgmesh->SubDomainsFromGeomSize[0]=Gh.NbSubDomains;
     760                bamgmesh->SubDomainsFromGeomSize[0]=Gh.nbsubdomains;
    762761                bamgmesh->SubDomainsFromGeomSize[1]=4;
    763                 if (Gh.NbSubDomains){
    764                         bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.NbSubDomains*sizeof(double));
    765                         for (i=0;i<Gh.NbSubDomains;i++){
     762                if (Gh.nbsubdomains){
     763                        bamgmesh->SubDomainsFromGeom=(double*)xmalloc(4*Gh.nbsubdomains*sizeof(double));
     764                        for (i=0;i<Gh.nbsubdomains;i++){
    766765                                bamgmesh->SubDomainsFromGeom[i*4+0]=2;
    767766                                bamgmesh->SubDomainsFromGeom[i*4+1]=Number(subdomains[i].edge)+1; //back to Matlab indexing
    768                                 bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].sens;
    769                                 bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ref;
     767                                bamgmesh->SubDomainsFromGeom[i*4+2]=subdomains[i].direction;
     768                                bamgmesh->SubDomainsFromGeom[i*4+3]=Gh.subdomains[i].ReferenceNumber;
    770769                        }
    771770                }
     
    806805                k=0;
    807806                for (i=0;i<nbe;i++){
    808                         if (edges[i].onGeometry) k=k+1;
     807                        if (edges[i].GeometricalEdgeHook) k=k+1;
    809808                }
    810809                bamgmesh->EdgesOnGeometricEdgeSize[0]=k;
     
    814813                        int count=0;
    815814                        for (i=0;i<nbe;i++){
    816                                 if (edges[i].onGeometry){
     815                                if (edges[i].GeometricalEdgeHook){
    817816                                        bamgmesh->EdgesOnGeometricEdge[count*2+0]=(double)i+1; //back to Matlab indexing
    818                                         bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].onGeometry)+1; //back to Matlab indexing
     817                                        bamgmesh->EdgesOnGeometricEdge[count*2+1]=(double)Gh.Number(edges[i].GeometricalEdgeHook)+1; //back to Matlab indexing
    819818                                        count=count+1;
    820819                                }
     
    14001399                                        edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]];
    14011400                                        edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]];
    1402                                         edges[add].onGeometry=NULL;
     1401                                        edges[add].GeometricalEdgeHook=NULL;
    14031402                                        //if already existed
    14041403                                        if (i<nbeold){
    1405                                                 edges[add].ref=edgessave[i].ref;                     
    1406                                                 edges[add].onGeometry=edgessave[i].onGeometry; //  HACK to get required edges
     1404                                                edges[add].ReferenceNumber=edgessave[i].ReferenceNumber;                     
     1405                                                edges[add].GeometricalEdgeHook=edgessave[i].GeometricalEdgeHook; //  HACK to get required edges
    14071406                                                printf("oh no...\n");
    14081407                                        }
    14091408                                        else
    1410                                          edges[add].ref=Min(edges[add].v[0]->ref(),edges[add].v[1]->ref());
     1409                                         edges[add].ReferenceNumber=Min(edges[add].v[0]->GetReferenceNumber(),edges[add].v[1]->GetReferenceNumber());
    14111410                                  }
    14121411                        }
     
    14741473                /*Reconstruct subdomains info*/
    14751474
    1476                 //check that NbSubDomains is empty
    1477                 if (NbSubDomains){
    1478                         ISSMERROR("NbSubDomains should be 0");
    1479                 }
    1480                 NbSubDomains=0;
     1475                //check that nbsubdomains is empty
     1476                if (nbsubdomains){
     1477                        ISSMERROR("nbsubdomains should be 0");
     1478                }
     1479                nbsubdomains=0;
    14811480
    14821481                //color the subdomains
     
    14941493
    14951494                                //color = number of subdomains
    1496                                 colorT[it]=NbSubDomains;
     1495                                colorT[it]=nbsubdomains;
    14971496
    14981497                                //color all the adjacent triangles of T that share a non marked edge
     
    15091508                                                //color the adjacent triangle
    15101509                                                if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) {
    1511                                                         colorT[jt]=NbSubDomains;
     1510                                                        colorT[jt]=nbsubdomains;
    15121511                                                        st[++level]=jt;
    15131512                                                        st[++level]=0;
     
    15171516                                        else level-=2;
    15181517                                }
    1519                                 NbSubDomains++;
    1520                         }
    1521                 }
    1522                 if (verbose> 3) printf("      The Number of sub domain = %i\n",NbSubDomains);
     1518                                nbsubdomains++;
     1519                        }
     1520                }
     1521                if (verbose> 3) printf("      The Number of sub domain = %i\n",nbsubdomains);
    15231522
    15241523                //build subdomains
    15251524                long isd;
    1526                 subdomains = new SubDomain[NbSubDomains];
     1525                subdomains = new SubDomain[nbsubdomains];
    15271526
    15281527                //initialize subdomains[isd].head as 0
    1529                 for (isd=0;isd<NbSubDomains;isd++) subdomains[isd].head =0;
     1528                for (isd=0;isd<nbsubdomains;isd++) subdomains[isd].head =0;
    15301529                 
    15311530                k=0;
     
    15351534                                if ((!tt || tt->color != triangles[it].color) && !subdomains[isd=colorT[it]].head){
    15361535                                        subdomains[isd].head = triangles+it;
    1537                                         subdomains[isd].ref =  triangles[it].color;
    1538                                         subdomains[isd].sens = j; // hack
     1536                                        subdomains[isd].ReferenceNumber =  triangles[it].color;
     1537                                        subdomains[isd].direction = j; // hack
    15391538                                        subdomains[isd].edge = 0;
    15401539                                        k++;
     
    15431542                }
    15441543                //check that we have been through all subdomains
    1545                 if (k!= NbSubDomains){
    1546                         ISSMERROR("k!= NbSubDomains");
     1544                if (k!= nbsubdomains){
     1545                        ISSMERROR("k!= nbsubdomains");
    15471546                }
    15481547                //delete colorT and st
     
    15691568                Gh.vertices = new GeometricalVertex[k];
    15701569                Gh.edges = new GeometricalEdge[nbe];
    1571                 Gh.NbSubDomains = NbSubDomains;
    1572                 Gh.subdomains = new GeometricalSubDomain[NbSubDomains];
     1570                Gh.nbsubdomains = nbsubdomains;
     1571                Gh.subdomains = new GeometricalSubDomain[nbsubdomains];
    15731572                if (verbose>3) printf("   number of vertices = %i\n   number of edges = %i\n",Gh.nbv,Gh.nbe);
    15741573                NbVerticesOnGeomVertex = Gh.nbv;
     
    16321631                        Gh.edges[i].tg[1]=R2();
    16331632
    1634                         bool required= edges[i].onGeometry;
     1633                        bool required= edges[i].GeometricalEdgeHook;
    16351634                        if(required) kreq++;
    1636                         edges[i].onGeometry =  Gh.edges + i;
     1635                        edges[i].GeometricalEdgeHook =  Gh.edges + i;
    16371636                        if(required){
    16381637                                Gh.edges[i].v[0]->SetRequired();
     
    16511650                        len[j1] += l12;
    16521651                        hmin = Min(hmin,l12);
    1653                         Gh.edges[i].ref  = edges[i].ref;
     1652                        Gh.edges[i].ReferenceNumber  = edges[i].ReferenceNumber;
    16541653
    16551654                        k = edge4->SortAndAdd(i0,i1);
     
    16701669
    16711670                //Build Gh.subdomains
    1672                 for (i=0;i<NbSubDomains;i++){
     1671                for (i=0;i<nbsubdomains;i++){
    16731672                        it = Number(subdomains[i].head);
    1674                         j = subdomains[i].sens;
     1673                        j = subdomains[i].direction;
    16751674                        long i0 = Number(triangles[it][VerticesOfTriangularEdge[j][0]]);
    16761675                        long i1 = Number(triangles[it][VerticesOfTriangularEdge[j][1]]);
    16771676                        k = edge4->SortAndFind(i0,i1);
    16781677                        if(k>=0){
    1679                                 subdomains[i].sens = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
     1678                                subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
    16801679                                subdomains[i].edge = edges+k;
    16811680                                Gh.subdomains[i].edge = Gh.edges + k;
    1682                                 Gh.subdomains[i].sens  =  subdomains[i].sens;
    1683                                 Gh.subdomains[i].ref =  subdomains[i].ref;
     1681                                Gh.subdomains[i].direction  =  subdomains[i].direction;
     1682                                Gh.subdomains[i].ReferenceNumber =  subdomains[i].ReferenceNumber;
    16841683                        }
    16851684                        else
     
    22032202                //  computed the number of cracked edge
    22042203                for (k=i=0;i<nbe;i++){
    2205                         if(edges[i].onGeometry->Cracked()) k++;
     2204                        if(edges[i].GeometricalEdgeHook->Cracked()) k++;
    22062205                }
    22072206
     
    22222221
    22232222                for (i=0;i<nbe;i++){
    2224                         if(edges[i].onGeometry->Cracked()){
     2223                        if(edges[i].GeometricalEdgeHook->Cracked()){
    22252224
    22262225                                //Fill edges fields of CrackedEdges
    2227                                 CrackedEdges[k  ].E =edges[i].onGeometry;
     2226                                CrackedEdges[k  ].E =edges[i].GeometricalEdgeHook;
    22282227                                CrackedEdges[k++].e1=&edges[i];
    22292228
     
    24592458
    24602459
    2461                         if (OutSide|| !Gh.subdomains || !Gh.NbSubDomains )
     2460                        if (OutSide|| !Gh.subdomains || !Gh.nbsubdomains )
    24622461                          { // No geom sub domain
    24632462                                long i;
    24642463                                if (subdomains) delete [] subdomains;
    24652464                                subdomains = new SubDomain[ NbSubDomTot];
    2466                                 NbSubDomains=  NbSubDomTot;
    2467                                 for ( i=0;i<NbSubDomains;i++) {
     2465                                nbsubdomains=  NbSubDomTot;
     2466                                for ( i=0;i<nbsubdomains;i++) {
    24682467                                        subdomains[i].head=NULL;
    2469                                         subdomains[i].ref=i+1;
     2468                                        subdomains[i].ReferenceNumber=i+1;
    24702469                                }
    24712470                                long * mark = new long[nbt];
     
    24882487                                                //    else if(mark[it] == -2 ) triangles[it].Draw(999);
    24892488                                                it++;} // end white (it<nbt)
    2490                                                 if (k!=NbSubDomains){
    2491                                                         ISSMERROR("k!=NbSubDomains");
     2489                                                if (k!=nbsubdomains){
     2490                                                        ISSMERROR("k!=nbsubdomains");
    24922491                                                }
    24932492                                                if(OutSide)
     
    24962495                                                        //  because in this case we have only the true boundary edge
    24972496                                                        // so teh boundary is manifold
    2498                                                         long nbk = NbSubDomains;
     2497                                                        long nbk = nbsubdomains;
    24992498                                                        while (nbk)
    25002499                                                         for (it=0;it<nbt && nbk ;it++)
     
    25052504                                                                  long kr = mark[it];
    25062505                                                                  if(kr !=kl) {
    2507                                                                           if (kl >=0 && subdomains[kl].ref <0 && kr >=0 && subdomains[kr].ref>=0)
    2508                                                                                 nbk--,subdomains[kr].ref=subdomains[kl].ref-1;
    2509                                                                           if (kr >=0 && subdomains[kr].ref <0 && kl >=0 && subdomains[kl].ref>=0)
    2510                                                                                 nbk--,subdomains[kl].ref=subdomains[kr].ref-1;
    2511                                                                           if(kr<0 && kl >=0 && subdomains[kl].ref>=0)
    2512                                                                                 nbk--,subdomains[kl].ref=-1;
    2513                                                                           if(kl<0 && kr >=0 && subdomains[kr].ref>=0)
    2514                                                                                 nbk--,subdomains[kr].ref=-1;
     2506                                                                          if (kl >=0 && subdomains[kl].ReferenceNumber <0 && kr >=0 && subdomains[kr].ReferenceNumber>=0)
     2507                                                                                nbk--,subdomains[kr].ReferenceNumber=subdomains[kl].ReferenceNumber-1;
     2508                                                                          if (kr >=0 && subdomains[kr].ReferenceNumber <0 && kl >=0 && subdomains[kl].ReferenceNumber>=0)
     2509                                                                                nbk--,subdomains[kl].ReferenceNumber=subdomains[kr].ReferenceNumber-1;
     2510                                                                          if(kr<0 && kl >=0 && subdomains[kl].ReferenceNumber>=0)
     2511                                                                                nbk--,subdomains[kl].ReferenceNumber=-1;
     2512                                                                          if(kl<0 && kr >=0 && subdomains[kr].ReferenceNumber>=0)
     2513                                                                                nbk--,subdomains[kr].ReferenceNumber=-1;
    25152514                                                                  }
    25162515                                                                 }
    25172516                                                        long  j=0;
    2518                                                         for ( i=0;i<NbSubDomains;i++)
    2519                                                          if((-subdomains[i].ref) %2) { // good
     2517                                                        for ( i=0;i<nbsubdomains;i++)
     2518                                                         if((-subdomains[i].ReferenceNumber) %2) { // good
    25202519                                                                 if(i != j)
    25212520                                                                  Exchange(subdomains[i],subdomains[j]);
     
    25302529                                                                 }//while (t)
    25312530                                                                }
    2532                                                         if(verbose>4) printf("      Number of removes subdomains (OutSideMesh) = %i\n",NbSubDomains-j);
    2533                                                         NbSubDomains=j;
     2531                                                        if(verbose>4) printf("      Number of removes subdomains (OutSideMesh) = %i\n",nbsubdomains-j);
     2532                                                        nbsubdomains=j;
    25342533                                                  }
    25352534
     
    25392538                        else
    25402539                          { // find the head for all sub domaine
    2541                                 if (Gh.NbSubDomains != NbSubDomains && subdomains)
     2540                                if (Gh.nbsubdomains != nbsubdomains && subdomains)
    25422541                                 delete [] subdomains, subdomains=0;
    25432542                                if (! subdomains  )
    2544                                  subdomains = new SubDomain[ Gh.NbSubDomains];
    2545                                 NbSubDomains =Gh.NbSubDomains;
     2543                                 subdomains = new SubDomain[ Gh.nbsubdomains];
     2544                                nbsubdomains =Gh.nbsubdomains;
    25462545                                long err=0;
    25472546                                CreateSingleVertexToTriangleConnectivity();
     
    25522551                                 mark[it]=triangles[it].link ? -1 : -2;
    25532552                                long inew =0;
    2554                                 for (int i=0;i<NbSubDomains;i++) {
     2553                                for (int i=0;i<nbsubdomains;i++) {
    25552554                                        GeometricalEdge &eg = *Gh.subdomains[i].edge;
    2556                                         subdomains[i].ref = Gh.subdomains[i].ref;
    2557                                         int ssdlab = subdomains[i].ref;
     2555                                        subdomains[i].ReferenceNumber = Gh.subdomains[i].ReferenceNumber;
     2556                                        int ssdlab = subdomains[i].ReferenceNumber;
    25582557                                        // by carefull is not easy to find a edge create from a GeometricalEdge
    25592558                                        // see routine MakeGeometricalEdgeToEdge
     
    25642563                                        BamgVertex * v0 =  e(0),*v1 = e(1);
    25652564                                        Triangle *t  = v0->t;
    2566                                         int sens = Gh.subdomains[i].sens;
    2567                                         // test if ge and e is in the same sens
    2568                                         if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) sens = -sens ;
    2569                                         subdomains[i].sens = sens;
     2565                                        int direction = Gh.subdomains[i].direction;
     2566                                        // test if ge and e is in the same direction
     2567                                        if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) direction = -direction ;
     2568                                        subdomains[i].direction = direction;
    25702569                                        subdomains[i].edge = &e;
    2571                                         if (!t || !sens){
    2572                                                 ISSMERROR("!t || !sens");
     2570                                        if (!t || !direction){
     2571                                                ISSMERROR("!t || !direction");
    25732572                                        }
    25742573
     
    25802579                                                }
    25812580                                                if (ta.EdgeVertex(0) == v1) { // ok we find the edge
    2582                                                         if (sens>0) 
     2581                                                        if (direction>0) 
    25832582                                                         subdomains[i].head=t=Adj(ta);
    25842583                                                        else
     
    26142613                                }
    26152614
    2616                                 if (inew < NbSubDomains) {
    2617                                         if (verbose>5) printf("WARNING: %i SubDomains are being removed\n",NbSubDomains-inew);
    2618                                         NbSubDomains=inew;}
     2615                                if (inew < nbsubdomains) {
     2616                                        if (verbose>5) printf("WARNING: %i SubDomains are being removed\n",nbsubdomains-inew);
     2617                                        nbsubdomains=inew;}
    26192618
    26202619
     
    26392638                srand(19999999);
    26402639                NbRef=0;
    2641                 NbOfTriangleSearchFind =0;
    2642                 NbOfSwapTriangle =0;
    26432640                nbv=0;
    26442641                maxnbv=maxnbv_in;
    26452642                nbt=0;
    2646                 NbOfQuad = 0;
     2643                nbq = 0;
    26472644                maxnbt=2*maxnbv_in-2;
    2648                 NbSubDomains=0;
     2645                nbsubdomains=0;
    26492646                NbVertexOnBThVertex=0;
    26502647                NbVertexOnBThEdge=0;
     
    26802677                NbVerticesOnGeomEdge=0;
    26812678                subdomains=NULL;
    2682                 NbSubDomains=0;
     2679                nbsubdomains=0;
    26832680        }
    26842681        /*}}}1*/
     
    29272924                  {
    29282925                        Edge * ei = edges+i;
    2929                         GeometricalEdge *onGeometry = ei->onGeometry;
    2930                         e[Gh.Number(onGeometry)] = ei;   
     2926                        GeometricalEdge *GeometricalEdgeHook = ei->GeometricalEdgeHook;
     2927                        e[Gh.Number(GeometricalEdgeHook)] = ei;   
    29312928                  }
    29322929                for ( i=0;i<nbe ; i++)
    29332930                 for (int ii=0;ii<2;ii++) {
    29342931                         Edge * ei = edges+i;
    2935                          GeometricalEdge *onGeometry = ei->onGeometry;
     2932                         GeometricalEdge *GeometricalEdgeHook = ei->GeometricalEdgeHook;
    29362933                         int j= ii;
    2937                          while (!(*onGeometry)[j].Required()) {
    2938                                  Adj(onGeometry,j); // next geom edge
     2934                         while (!(*GeometricalEdgeHook)[j].Required()) {
     2935                                 Adj(GeometricalEdgeHook,j); // next geom edge
    29392936                                 j=1-j;
    2940                                  if (e[Gh.Number(onGeometry)])  break; // optimisation
    2941                                  e[Gh.Number(onGeometry)] = ei;
     2937                                 if (e[Gh.Number(GeometricalEdgeHook)])  break; // optimisation
     2938                                 e[Gh.Number(GeometricalEdgeHook)] = ei;
    29422939                         }
    29432940                 }
     
    29892986                                 triangles[i].SetHidden(j),kk++;
    29902987                          }
    2991                         NbOfQuad = kk;
     2988                        nbq = kk;
    29922989                        if (verbose>2){
    2993                                 printf("   number of quadrilaterals    = %i\n",NbOfQuad);
    2994                                 printf("   number of triangles         = %i\n",nbt-nbtout- NbOfQuad*2);
     2990                                printf("   number of quadrilaterals    = %i\n",nbq);
     2991                                printf("   number of triangles         = %i\n",nbt-nbtout- nbq*2);
    29952992                                printf("   number of outside triangles = %i\n",nbtout);
    29962993                        }
     
    31093106                        for (i=0;i<Bh.nbv;i++){
    31103107                                BamgVertex &bv=Bh[i];
    3111                                 if (!bv.onGeometry){
     3108                                if (!bv.GeometricalEdgeHook){
    31123109                                        vertices[nbv].r   = bv.r;
    31133110                                        vertices[nbv++].m = bv.m;
     
    32183215                BamgVertex * pvA=&vA, * pvB=&vB;
    32193216                if (vA.vint == IsVertexOnVertex){
    3220                         pA=vA.onBackgroundVertex;
     3217                        pA=vA.BackgroundVertexHook;
    32213218                }
    32223219                else if (vA.vint == IsVertexOnEdge){
    3223                         pA=vA.onBackgroundEdge->be;
    3224                         tA=vA.onBackgroundEdge->abcisse;
     3220                        pA=vA.BackgroundEdgeHook->be;
     3221                        tA=vA.BackgroundEdgeHook->abcisse;
    32253222                }
    32263223                else {
     
    32293226
    32303227                if (vB.vint == IsVertexOnVertex){
    3231                         pB=vB.onBackgroundVertex;
     3228                        pB=vB.BackgroundVertexHook;
    32323229                }
    32333230                else if(vB.vint == IsVertexOnEdge){
    3234                         pB=vB.onBackgroundEdge->be;
    3235                         tB=vB.onBackgroundEdge->abcisse;
     3231                        pB=vB.BackgroundEdgeHook->be;
     3232                        tB=vB.BackgroundEdgeHook->abcisse;
    32363233                }
    32373234                else {
     
    32553252                if( vA.vint == IsVertexOnEdge)
    32563253                  { // find the start edge
    3257                         e = vA.onBackgroundEdge->be;     
     3254                        e = vA.BackgroundEdgeHook->be;   
    32583255
    32593256                  }
     
    32653262                        Exchange(pvA,pvB);
    32663263                        Exchange(A,B);
    3267                         e =  vB.onBackgroundEdge->be;
     3264                        e =  vB.BackgroundEdgeHook->be;
    32683265
    32693266                  }
     
    32723269                  }
    32733270
    3274                 // find the direction of walking with sens of edge and pA,PB;
     3271                // find the direction of walking with direction of edge and pA,PB;
    32753272                R2 AB=B-A;
    32763273
    32773274                double cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB);
    32783275                int kkk=0;
    3279                 int sens = (cosE01AB>0) ? 1 : 0;
     3276                int direction = (cosE01AB>0) ? 1 : 0;
    32803277
    32813278                //   double l=0; // length of the edge AB
     
    32923289                        double te0;
    32933290                        // we suppose take the curve's abcisse
    3294                         for ( eee=e,iii=sens,te0=tA;
     3291                        for ( eee=e,iii=direction,te0=tA;
    32953292                                                eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ;
    32963293                                                neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) {
     
    33623359        for (i=0;i<nbv;i++) ordre[i]=0;
    33633360
    3364         //Initialize NbSubDomains
    3365         NbSubDomains =0;
     3361        //Initialize nbsubdomains
     3362        nbsubdomains =0;
    33663363
    33673364        /* generation of triangles adjacency*/
     
    36193616        for (i=0;i<nbe;i++){
    36203617        /*If the current mesh edge is on Geometry*/
    3621                 if(edges[i].onGeometry){
     3618                if(edges[i].GeometricalEdgeHook){
    36223619                        for(int j=0;j<2;j++){
    36233620                                /*Go through the edges adjacent to current edge (if on the same curve)*/
     
    36253622                                        /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/
    36263623                                        /*Check that the 2 vertices are on geometry AND required*/
    3627                                         if(!edges[i][j].onGeometry->IsRequiredVertex()){
     3624                                        if(!edges[i][j].GeometricalEdgeHook->IsRequiredVertex()){
    36283625                                                printf("ReconstructExistingMesh error message: problem with the edge number %i: [%i %i]\n",i+1,Number(edges[i][0])+1,Number(edges[i][1])+1);
    3629                                                 printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].onGeometry)+1);
    3630                                                 if (edges[i][j].onGeometry->OnGeomVertex())
    3631                                                  printf("the vertex number %i of this edge is a geometric BamgVertex number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->gv)+1);
    3632                                                 else if (edges[i][j].onGeometry->OnGeomEdge())
    3633                                                  printf("the vertex number %i of this edge is a geometric Edge number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].onGeometry->ge)+1);
     3626                                                printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].GeometricalEdgeHook)+1);
     3627                                                if (edges[i][j].GeometricalEdgeHook->OnGeomVertex())
     3628                                                 printf("the vertex number %i of this edge is a geometric BamgVertex number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].GeometricalEdgeHook->gv)+1);
     3629                                                else if (edges[i][j].GeometricalEdgeHook->OnGeomEdge())
     3630                                                 printf("the vertex number %i of this edge is a geometric Edge number %i\n",Number(edges[i][j])+1,Gh.Number(edges[i][j].GeometricalEdgeHook->ge)+1);
    36343631                                                else
    3635                                                  printf("Its pointer is %p\n",edges[i][j].onGeometry);
     3632                                                 printf("Its pointer is %p\n",edges[i][j].GeometricalEdgeHook);
    36363633
    36373634                                                printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");
     
    36553652                for ( it=0;it<nbt;it++)
    36563653                 renu[it]=-1; // outside triangle
    3657                 for ( i=0;i<NbSubDomains;i++)
     3654                for ( i=0;i<nbsubdomains;i++)
    36583655                  {
    36593656                        t=t0=subdomains[i].head;
     
    36903687                 triangles[it].Renumbering(triangles,te,renu);
    36913688
    3692                 for ( i=0;i<NbSubDomains;i++)
     3689                for ( i=0;i<nbsubdomains;i++)
    36933690                 subdomains[i].head=triangles+renu[Number(subdomains[i].head)];
    36943691
     
    39893986                 if (tstart[i] != &vide) // not a boundary vertex
    39903987                  delta=Max(delta,vertices[i].Smoothing(*this,BTh,tstart[i],omega));
    3991                 if (!NbOfQuad)
     3988                if (!nbq)
    39923989                 for ( i=0;i<nbv;i++)
    39933990                  if (tstart[i] != &vide) // not a boundary vertex
     
    41004097                long newnbt=0,newnbv=0;
    41014098                long * kedge = 0;
    4102                 long newNbOfQuad=NbOfQuad;
     4099                long newnbq=nbq;
    41034100                long * ksplit= 0, * ksplitarray=0;
    41044101                long kkk=0;
     
    41174114                Triangle * lastT = triangles + nbt;
    41184115                for (i=0;i<nbe;i++)
    4119                  if(edges[i].onGeometry) NbEdgeOnGeom++;
     4116                 if(edges[i].GeometricalEdgeHook) NbEdgeOnGeom++;
    41204117                long newnbe=nbe+nbe;
    41214118                //  long newNbVerticesOnGeomVertex=NbVerticesOnGeomVertex;
     
    41444141                long ferr=0;
    41454142                for (i=0;i<nbe;i++)
    4146                  newedges[ie].onGeometry=0;
     4143                 newedges[ie].GeometricalEdgeHook=0;
    41474144
    41484145                for (i=0;i<nbe;i++)
    41494146                  {
    4150                         GeometricalEdge *ong =  edges[i].onGeometry;
     4147                        GeometricalEdge *ong =  edges[i].GeometricalEdgeHook;
    41514148
    41524149                        newedges[ie]=edges[i];
     
    41674164                                                ISSMERROR("!edgesGtoB");
    41684165                                        }
    4169                                         ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].onGeometry)],
     4166                                        ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].GeometricalEdgeHook)],
    41704167                                                                edges[i][0],edges[i][1],0.5,vertices[k],
    41714168                                                                newVertexOnBThEdge[kvb],
    41724169                                                                newVerticesOnGeomEdge[kvg++]);
    4173                                         vertices[k].ReferenceNumber= edges[i].ref;
     4170                                        vertices[k].ReferenceNumber= edges[i].ReferenceNumber;
    41744171                                        vertices[k].DirOfSearch =   NoDirOfSearch;       
    41754172                                        ;
     
    41874184                                                                0.5,vertices[k],newVerticesOnGeomEdge[kvg++]);
    41884185                                        // vertices[k].i = toI2( vertices[k].r);
    4189                                         vertices[k].ReferenceNumber = edges[i].ref;
     4186                                        vertices[k].ReferenceNumber = edges[i].ReferenceNumber;
    41904187                                        vertices[k].DirOfSearch = NoDirOfSearch;
    41914188                                        vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);             
     
    41964193                                vertices[k].r = ((R2) edges[i][0] + (R2)  edges[i][1] )*0.5;
    41974194                                vertices[k].m =  Metric(0.5,edges[i][0],0.5,edges[i][1]);
    4198                                 vertices[k].onGeometry = 0;
     4195                                vertices[k].GeometricalEdgeHook = 0;
    41994196                          }
    42004197                        //vertices[k].i = toI2( vertices[k].r);
     
    42024199                        R2 AA = (A+AB)*0.5;
    42034200                        R2 BB = (AB+B)*0.5;
    4204                         vertices[k].ReferenceNumber = edges[i].ref;
     4201                        vertices[k].ReferenceNumber = edges[i].ReferenceNumber;
    42054202                        vertices[k].DirOfSearch = NoDirOfSearch;
    42064203
    4207                         newedges[ie].onGeometry = Gh.Containing(AA,ong);
     4204                        newedges[ie].GeometricalEdgeHook = Gh.Containing(AA,ong);
    42084205                        newedges[ie++].v[1]=vertices+k;
    42094206
     
    42114208                        newedges[ie].adj[0]=newedges + ie -1;
    42124209                        newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ;
    4213                         newedges[ie].onGeometry =  Gh.Containing(BB,ong);
     4210                        newedges[ie].GeometricalEdgeHook =  Gh.Containing(BB,ong);
    42144211                        newedges[ie++].v[0]=vertices+k;
    42154212                        k++;
     
    43494346                  }
    43504347                //  now do the element split
    4351                 newNbOfQuad = 4*NbOfQuad;
     4348                newnbq = 4*nbq;
    43524349                nbv = k;
    43534350                kkk = nbt;
     
    45174514                        }
    45184515
    4519                         if (kk==6)  newNbOfQuad+=3;
     4516                        if (kk==6)  newnbq+=3;
    45204517                        for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk;
    45214518                        ksplit[i]= nbt; // save last adresse of the new triangles
     
    45444541                edges = newedges;
    45454542                nbe = newnbe;
    4546                 NbOfQuad = newNbOfQuad;
    4547 
    4548                 for (i=0;i<NbSubDomains;i++)
     4543                nbq = newnbq;
     4544
     4545                for (i=0;i<nbsubdomains;i++)
    45494546                  {
    45504547                        long k = subdomains[i].edge- edges;
     
    45664563
    45674564                if (verbose>2){
    4568                         printf("   number of quadrilaterals    = %i\n",NbOfQuad);
    4569                         printf("   number of triangles         = %i\n",nbt-nbtout- NbOfQuad*2);
     4565                        printf("   number of quadrilaterals    = %i\n",nbq);
     4566                        printf("   number of triangles         = %i\n",nbt-nbtout- nbq*2);
    45704567                        printf("   number of outside triangles = %i\n",nbtout);
    45714568                }
     
    46074604                                  BamgVertex &v0 = t[VerticesOfTriangularEdge[j][0]];
    46084605                                  BamgVertex &v1 = t[VerticesOfTriangularEdge[j][1]];
    4609                                   if (v0.onGeometry && v1.onGeometry){
     4606                                  if (v0.GeometricalEdgeHook && v1.GeometricalEdgeHook){
    46104607                                          R2 P= ((R2) v0 + (R2) v1)*0.5;
    46114608                                          if ( nbv<maxnbv) {
     
    47504747        if (t->det<0) // outside triangle
    47514748         dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;
    4752         //  NbOfTriangleSearchFind += counter; 
    47534749        return t;
    47544750}
     
    47774773
    47784774        //loop over all subdomains
    4779         for (int i=0;i<NbSubDomains;i++){
     4775        for (int i=0;i<nbsubdomains;i++){
    47804776
    47814777                //first triangle of the subdomain i
     
    48144810
    48154811                int i,j,k;
    4816                 int NbOfCurves=0,NbNewPoints,NbEdgeCurve;
     4812                int nbcurves=0,NbNewPoints,NbEdgeCurve;
    48174813                double lcurve,lstep,s;
    48184814                const int MaxSubEdge = 10;
     
    48904886                        long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
    48914887                        Gh.UnMarkEdges();       
    4892                         NbOfCurves=0;
     4888                        nbcurves=0;
    48934889
    48944890                        //go through the edges of the geometry
     
    49224918                                                                        edges[nbe].v[0]=a->to;
    49234919                                                                        edges[nbe].v[1]=b->to;;
    4924                                                                         edges[nbe].ref = e->ref;
    4925                                                                         edges[nbe].onGeometry = e;
     4920                                                                        edges[nbe].ReferenceNumber = e->ReferenceNumber;
     4921                                                                        edges[nbe].GeometricalEdgeHook = e;
    49264922                                                                        edges[nbe].adj[0] = 0;
    49274923                                                                        edges[nbe].adj[1] = 0;
     
    50215017                                                                                vb = &vertices[nbv++];
    50225018                                                                                vb->m = Metric(aa,a->m,bb,b->m);
    5023                                                                                 vb->ReferenceNumber = e->ref;
     5019                                                                                vb->ReferenceNumber = e->ReferenceNumber;
    50245020                                                                                vb->DirOfSearch =NoDirOfSearch;
    50255021                                                                                double abcisse = k ? bb : aa;
     
    50315027                                                                                edges[nbe].v[0]=va;
    50325028                                                                                edges[nbe].v[1]=vb;
    5033                                                                                 edges[nbe].ref =e->ref;
    5034                                                                                 edges[nbe].onGeometry = e;
     5029                                                                                edges[nbe].ReferenceNumber =e->ReferenceNumber;
     5030                                                                                edges[nbe].GeometricalEdgeHook = e;
    50355031                                                                                edges[nbe].adj[0] = PreviousNewEdge;
    50365032                                                                                if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
     
    50535049                                                                if(!kstep){
    50545050                                                                        NbVerticesOnGeomEdge0 += NbNewPoints;
    5055                                                                         NbOfCurves++;
     5051                                                                        nbcurves++;
    50565052                                                                }
    50575053                                                                nbvend=nbv+NbNewPoints;
     
    50615057                                                                edges[nbe].v[0]=va;
    50625058                                                                edges[nbe].v[1]=vb;
    5063                                                                 edges[nbe].ref = e->ref;
    5064                                                                 edges[nbe].onGeometry = e;
     5059                                                                edges[nbe].ReferenceNumber = e->ReferenceNumber;
     5060                                                                edges[nbe].GeometricalEdgeHook = e;
    50655061                                                                edges[nbe].adj[0] = PreviousNewEdge;
    50665062                                                                edges[nbe].adj[1] = 0;
     
    51465142                this->Init(imaxnbv);
    51475143                BTh.SetVertexFieldOn();
    5148                 int* bcurve = new int[Gh.NbOfCurves]; //
     5144                int* bcurve = new int[Gh.nbcurves]; //
    51495145
    51505146                /* There are 2 ways to make the loop
     
    51965192                Gh.UnMarkEdges();       
    51975193                int bfind=0;
    5198                 for (int i=0;i<Gh.NbOfCurves;i++) bcurve[i]=-1;
     5194                for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1;
    51995195
    52005196                /*Loop over the backgrounf mesh BTh edges*/
     
    52065202
    52075203                                /* If one of the vertex is required we are in a new curve*/
    5208                                 if (ei[je].onGeometry->IsRequiredVertex()){
     5204                                if (ei[je].GeometricalEdgeHook->IsRequiredVertex()){
    52095205
    52105206                                        /*Get curve number*/
    5211                                         int nc=ei.onGeometry->CurveNumber;
     5207                                        int nc=ei.GeometricalEdgeHook->CurveNumber;
    52125208                                       
    52135209                                        //printf("Dealing with curve number %i\n",nc);
    5214                                         //printf("edge on geometry is same as GhCurve? %s\n",(ei.onGeometry==Gh.curves[nc].be || ei.onGeometry==Gh.curves[nc].ee)?"yes":"no");
    5215                                         //if(ei.onGeometry==Gh.curves[nc].be || ei.onGeometry==Gh.curves[nc].ee){
    5216                                         //      printf("Do we have the right extremity? curve first vertex -> %s\n",((GeometricalVertex *)*ei[je].onGeometry==&(*Gh.curves[nc].be)[Gh.curves[nc].kb])?"yes":"no");
    5217                                         //      printf("Do we have the right extremity? curve last  vertex -> %s\n",((GeometricalVertex *)*ei[je].onGeometry==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke])?"yes":"no");
     5210                                        //printf("edge on geometry is same as GhCurve? %s\n",(ei.GeometricalEdgeHook==Gh.curves[nc].be || ei.GeometricalEdgeHook==Gh.curves[nc].ee)?"yes":"no");
     5211                                        //if(ei.GeometricalEdgeHook==Gh.curves[nc].be || ei.GeometricalEdgeHook==Gh.curves[nc].ee){
     5212                                        //      printf("Do we have the right extremity? curve first vertex -> %s\n",((GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].be)[Gh.curves[nc].kb])?"yes":"no");
     5213                                        //      printf("Do we have the right extremity? curve last  vertex -> %s\n",((GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke])?"yes":"no");
    52185214                                        //}
    52195215                                        //BUG FIX from original bamg
    52205216                                        /*Check that we are on the same edge and right vertex (0 or 1) */
    5221                                         if(ei.onGeometry==Gh.curves[nc].be  && (GeometricalVertex *)*ei[je].onGeometry==&(*Gh.curves[nc].be)[Gh.curves[nc].kb]){
     5217                                        if(ei.GeometricalEdgeHook==Gh.curves[nc].be  && (GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].be)[Gh.curves[nc].kb]){
    52225218                                                bcurve[nc]=iedge*2+je;
    52235219                                                bfind++;       
    52245220                                        }
    5225                                         else if ((ei.onGeometry==Gh.curves[nc].ee  && (GeometricalVertex *)*ei[je].onGeometry==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke]) && bcurve[nc]==-1){
     5221                                        else if ((ei.GeometricalEdgeHook==Gh.curves[nc].ee  && (GeometricalVertex *)*ei[je].GeometricalEdgeHook==&(*Gh.curves[nc].ee)[Gh.curves[nc].ke]) && bcurve[nc]==-1){
    52265222                                                bcurve[nc]=iedge*2+je;
    52275223                                                bfind++;       
     
    52305226                        }
    52315227                }
    5232                 if (bfind!=Gh.NbOfCurves) ISSMERROR("problem generating number of curves (Gh.NbOfCurves=%i bfind=%i)",Gh.NbOfCurves,bfind);
     5228                if (bfind!=Gh.nbcurves) ISSMERROR("problem generating number of curves (Gh.nbcurves=%i bfind=%i)",Gh.nbcurves,bfind);
    52335229
    52345230                // method in 2 + 1 step
     
    52475243
    52485244                        /*Go through all geometrical curve*/
    5249                         for (int icurve=0;icurve<Gh.NbOfCurves;icurve++){
     5245                        for (int icurve=0;icurve<Gh.nbcurves;icurve++){
    52505246
    52515247                                /*Get edge and vertex (index) of background mesh on this curve*/
     
    52825278                                                int k0equi=jedgeequi,k1equi;             
    52835279                                                Edge * peequi= BTh.edges+iedgeequi;
    5284                                                 GeometricalEdge *ongequi = peequi->onGeometry;
     5280                                                GeometricalEdge *ongequi = peequi->GeometricalEdgeHook;
    52855281
    52865282                                                double sNew=Lstep;// abscisse of the new points (phase==1)
    52875283                                                L=0;// length of the curve
    52885284                                                long i=0;// index of new points on the curve
    5289                                                 register GeometricalVertex * GA0 = *(*peequi)[k0equi].onGeometry;
     5285                                                register GeometricalVertex * GA0 = *(*peequi)[k0equi].GeometricalEdgeHook;
    52905286                                                BamgVertex *A0;
    52915287                                                A0 = GA0->to;  // the vertex in new mesh
     
    52975293                                                ISSMASSERT(A0-vertices>=0 && A0-vertices<nbv);
    52985294                                                if(ongequi->Required()){
    5299                                                         GeometricalVertex *GA1 = *(*peequi)[1-k0equi].onGeometry;
     5295                                                        GeometricalVertex *GA1 = *(*peequi)[1-k0equi].GeometricalEdgeHook;
    53005296                                                        A1 = GA1->to;  //
    53015297                                                }       
     
    53065302                                                                k1 = 1-k0; // next vertex of the edge
    53075303                                                                k1equi= 1 - k0equi;
    5308                                                                 ISSMASSERT(pe && ee.onGeometry);
    5309                                                                 ee.onGeometry->SetMark();
     5304                                                                ISSMASSERT(pe && ee.GeometricalEdgeHook);
     5305                                                                ee.GeometricalEdgeHook->SetMark();
    53105306                                                                BamgVertex & v0=ee[0], & v1=ee[1];
    53115307                                                                R2 AB=(R2)v1-(R2)v0;
     
    53415337                                                                                VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
    53425338                                                                                ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
    5343                                                                                 A1->ReferenceNumber = eeequi.ref;
     5339                                                                                A1->ReferenceNumber = eeequi.ReferenceNumber;
    53445340                                                                                A1->DirOfSearch =NoDirOfSearch;
    5345                                                                                 e->onGeometry = ongequi;
     5341                                                                                e->GeometricalEdgeHook = ongequi;
    53465342                                                                                e->v[0]=A0;
    53475343                                                                                e->v[1]=A1;
    5348                                                                                 e->ref = eeequi.ref;
     5344                                                                                e->ReferenceNumber = eeequi.ReferenceNumber;
    53495345                                                                                e->adj[0]=PreviousNewEdge;
    53505346
     
    53585354
    53595355                                                                //some checks
    5360                                                                 ISSMASSERT(ee.onGeometry->CurveNumber==ei.onGeometry->CurveNumber);
    5361                                                                 if (ee[k1].onGeometry->IsRequiredVertex()) {
    5362                                                                         ISSMASSERT(eeequi[k1equi].onGeometry->IsRequiredVertex());
    5363                                                                         register GeometricalVertex * GA1 = *eeequi[k1equi].onGeometry;
     5356                                                                ISSMASSERT(ee.GeometricalEdgeHook->CurveNumber==ei.GeometricalEdgeHook->CurveNumber);
     5357                                                                if (ee[k1].GeometricalEdgeHook->IsRequiredVertex()) {
     5358                                                                        ISSMASSERT(eeequi[k1equi].GeometricalEdgeHook->IsRequiredVertex());
     5359                                                                        register GeometricalVertex * GA1 = *eeequi[k1equi].GeometricalEdgeHook;
    53645360                                                                        A1=GA1->to;// the vertex in new mesh
    53655361                                                                        ISSMASSERT(A1-vertices>=0 && A1-vertices<nbv);
     
    53795375                                                if (phase){ // construction of the last edge
    53805376                                                        Edge* e=edges + nbe++;
    5381                                                         e->onGeometry  = ongequi;
     5377                                                        e->GeometricalEdgeHook  = ongequi;
    53825378                                                        e->v[0]=A0;
    53835379                                                        e->v[1]=A1;
    5384                                                         e->ref = peequi->ref;
     5380                                                        e->ReferenceNumber = peequi->ReferenceNumber;
    53855381                                                        e->adj[0]=PreviousNewEdge;
    53865382                                                        e->adj[1]=0;
     
    56495645        AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.vint][0]);
    56505646        BamgVertex   *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
    5651         // we turn around a in the  direct sens 
     5647        // we turn around a in the  direct direction 
    56525648
    56535649        Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1;
     
    58165812                        if(!ToSwap) tt1 =  Next(tt2);
    58175813                }
    5818                 else { // changement de sens
     5814                else { // changement de direction
    58195815                        ret = -1;
    58205816                        Exchange(pva,pvb);
  • issm/trunk/src/c/objects/Bamg/Mesh.h

    r5146 r5148  
    3434                        SubDomain                    *subdomains;
    3535                        long                          NbRef;                 // counter of ref on the this class if 0 we can delete
    36                         long                          maxnbv,maxnbt;           // nombre max de sommets , de triangles
    37                         long                          nbv,nbt,nbe;           // nb of vertices, of triangles, of edges
    38                         long                          NbOfQuad;              // nb of quadrangle
    39                         long                          NbSubDomains;
     36                        long                          maxnbv,maxnbt;         // nombre max de sommets , de triangles
     37                        long                          nbv,nbt,nbe,nbq;       // nb of vertices, of triangles, of edges and quadrilaterals
     38                        long                          nbsubdomains;
    4039                        long                          nbtout;                // Nb of oudeside triangle
    41                         long                          NbOfTriangleSearchFind;
    42                         long                          NbOfSwapTriangle;
     40
     41                        R2                            pmin,pmax;             // extrema
     42                        double                        coefIcoor;             // coef to integer Icoor1;
     43                        ListofIntersectionTriangles   lIntTria;
     44
    4345                        long                          NbVerticesOnGeomVertex;
    4446                        VertexOnGeom                 *VerticesOnGeomVertex;
     
    5355                        long                          NbCrackedEdges;
    5456                        CrackedEdge                  *CrackedEdges;
    55                         R2                            pmin,pmax;             // extrema
    56                         double                        coefIcoor;             // coef to integer Icoor1;
    57                         ListofIntersectionTriangles   lIntTria;
    5857
    5958                        //Constructors/Destructors
     
    139138                          }
    140139                        inline  void  SetVertexFieldOn(){
    141                                 for (int i=0;i<nbv;i++)                    vertices[i].onGeometry=NULL;
     140                                for (int i=0;i<nbv;i++)                    vertices[i].GeometricalEdgeHook=NULL;
    142141                                for (int j=0;j<NbVerticesOnGeomVertex;j++) VerticesOnGeomVertex[j].SetOn();
    143142                                for (int k=0;k<NbVerticesOnGeomEdge;k++ )  VerticesOnGeomEdge[k].SetOn();
    144143                        }             
    145144                        inline  void   SetVertexFieldOnBTh(){
    146                                 for (int i=0;i<nbv;i++)                 vertices[i].onGeometry=NULL;
     145                                for (int i=0;i<nbv;i++)                 vertices[i].GeometricalEdgeHook=NULL;
    147146                                for (int j=0;j<NbVertexOnBThVertex;j++) VertexOnBThVertex[j].SetOnBTh();
    148147                                for (int k=0;k<NbVertexOnBThEdge;k++ )  VertexOnBThEdge[k].SetOnBTh();
  • issm/trunk/src/c/objects/Bamg/Metric.cpp

    r5016 r5148  
    156156                                  {
    157157                                        k=k/2;
    158                                         // we begin by the end to walk in the correct sens from a to b
     158                                        // we begin by the end to walk in the correct direction from a to b
    159159                                        // due to the stack
    160160                                        Ms1[level]=Mi;
  • issm/trunk/src/c/objects/Bamg/QuadTree.cpp

    r5120 r5148  
    388388                          {
    389389                                I2 i2 =  b->v[k]->i;
    390                                 //   try if is in the right sens --
     390                                //   try if is in the right direction --
    391391                                h0 = NORM(iplus,i2.x,jplus,i2.y);
    392392                                if (h0 <h) {
     
    415415                                        NbVerticesSearch++;
    416416                                        I2 i2 =  b->v[k]->i;
    417                                         // if good sens when try --
     417                                        // if good direction when try --
    418418
    419419                                        h0 = NORM(iplus,i2.x,jplus,i2.y);
  • issm/trunk/src/c/objects/Bamg/R2.h

    r5124 r5148  
    9191          }
    9292        template  <class R,class RR> 
    93           inline R NormeInfini (const P2<R,RR> x) {
    94                   return Max(Abs(x.x),Abs(x.y)) ;
    95           }
    96         template  <class R,class RR> 
    9793          inline RR Norme2_2 (const P2<R,RR> x) {
    9894                  return (RR)x.x*(RR)x.x + (RR)x.y*(RR)x.y ;
  • issm/trunk/src/c/objects/Bamg/SubDomain.h

    r5095 r5148  
    1212
    1313        class SubDomain {
     14
    1415                public:
    15                         Triangle * head;
    16                         long  ref; 
    17                         int sens; // -1 or 1
    18                         Edge* edge; // to  geometrical 
     16
     17                        Triangle *head;
     18                        long      ReferenceNumber;
     19                        int       direction;   // -1 or 1
     20                        Edge     *edge;        // to geometrical
    1921
    2022                        //Methods
  • issm/trunk/src/c/objects/Bamg/Triangle.cpp

    r5143 r5148  
    1919                        ISSMERROR("i>=nbv || j>=nbv || k>=nbv");
    2020                }
    21                 TriaVertices[0]=v+i;
    22                 TriaVertices[1]=v+j;
    23                 TriaVertices[2]=v+k;
    24                 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
    25                 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
     21                vertices[0]=v+i;
     22                vertices[1]=v+j;
     23                vertices[2]=v+k;
     24                adj[0]=adj[1]=adj[2]=0;
     25                AdjEdgeNumber[0]=AdjEdgeNumber[1]=AdjEdgeNumber[2]=0;
    2626                det=0;
    2727        }
     
    2929        /*FUNCTION Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2) {{{1*/
    3030        Triangle::Triangle(BamgVertex *v0,BamgVertex *v1,BamgVertex *v2){
    31                 TriaVertices[0]=v0;
    32                 TriaVertices[1]=v1;
    33                 TriaVertices[2]=v2;
    34                 TriaAdjTriangles[0]=TriaAdjTriangles[1]=TriaAdjTriangles[2]=0;
    35                 TriaAdjSharedEdge[0]=TriaAdjSharedEdge[1]=TriaAdjSharedEdge[2]=0;
     31                vertices[0]=v0;
     32                vertices[1]=v1;
     33                vertices[2]=v2;
     34                adj[0]=adj[1]=adj[2]=0;
     35                AdjEdgeNumber[0]=AdjEdgeNumber[1]=AdjEdgeNumber[2]=0;
    3636                if (v0) det=0;
    3737                else {
     
    4848
    4949                printf("Triangle:\n");
    50                 printf("   TriaVertices pointer towards three vertices\n");
    51                 printf("      TriaVertices[0] TriaVertices[1] TriaVertices[2] = %p %p %p\n",TriaVertices[0],TriaVertices[1],TriaVertices[2]);
    52                 printf("   TriaAdjTriangles pointer towards three adjacent triangles\n");
    53                 printf("      TriaAdjTriangles[0] TriaAdjTriangles[1] TriaAdjTriangles[2] = %p %p %p\n",TriaAdjTriangles[0],TriaAdjTriangles[1],TriaAdjTriangles[2]);
     50                printf("   vertices pointer towards three vertices\n");
     51                printf("      vertices[0] vertices[1] vertices[2] = %p %p %p\n",vertices[0],vertices[1],vertices[2]);
     52                printf("   adj pointer towards three adjacent triangles\n");
     53                printf("      adj[0] adj[1] adj[2] = %p %p %p\n",adj[0],adj[1],adj[2]);
    5454                printf("   det (integer triangle determinant) = %i\n",det);
    5555                if (link){
     
    6262                printf("\nThree vertices:\n");
    6363                for(i=0;i<3;i++){
    64                         if (TriaVertices[i]){
    65                                 TriaVertices[i]->Echo();
     64                        if (vertices[i]){
     65                                vertices[i]->Echo();
    6666                        }
    6767                        else{
     
    9999                        k++;
    100100                        //Get ttc, adjacent triangle of t with respect to vertex j
    101                         ttc =  t->TriaAdjTriangles[j];
     101                        ttc =  t->adj[j];
    102102                        //is the current triangle inside or outside?
    103103                        outside = !ttc->link;
     
    108108                        t = ttc;
    109109                        //NextEdge[3] = {1,2,0};
    110                         jc = NextEdge[t->TriaAdjSharedEdge[j]&3];
     110                        jc = NextEdge[t->AdjEdgeNumber[j]&3];
    111111                        j = NextEdge[jc];
    112112
     
    132132
    133133                // initialize tp, jp the previous triangle & edge
    134                 Triangle *tp=TriaAdjTriangles[jp];
    135                 jp = TriaAdjSharedEdge[jp]&3;
     134                Triangle *tp=adj[jp];
     135                jp = AdjEdgeNumber[jp]&3;
    136136                do {
    137137                        while (t->swap(j,koption)){
     
    139139                                NbSwap++;
    140140                                k++;
    141                                 t=  tp->TriaAdjTriangles[jp];      // set unchange t qnd j for previous triangles
    142                                 j=  NextEdge[tp->TriaAdjSharedEdge[jp]&3];
     141                                t=  tp->adj[jp];      // set unchange t qnd j for previous triangles
     142                                j=  NextEdge[tp->AdjEdgeNumber[jp]&3];
    143143                        }
    144144                        // end on this  Triangle
     
    146146                        jp = NextEdge[j];
    147147
    148                         t=  tp->TriaAdjTriangles[jp];      // set unchange t qnd j for previous triangles
    149                         j=  NextEdge[tp->TriaAdjSharedEdge[jp]&3];
     148                        t=  tp->adj[jp];      // set unchange t qnd j for previous triangles
     149                        j=  NextEdge[tp->AdjEdgeNumber[jp]&3];
    150150
    151151                } while( t != this);
     
    159159                if (link) {
    160160                        int a=-1;
    161                         if (TriaAdjSharedEdge[0] & 16 ) a=0;
    162                         if (TriaAdjSharedEdge[1] & 16 ) a=1;
    163                         if (TriaAdjSharedEdge[2] & 16 ) a=2;
     161                        if (AdjEdgeNumber[0] & 16 ) a=0;
     162                        if (AdjEdgeNumber[1] & 16 ) a=1;
     163                        if (AdjEdgeNumber[2] & 16 ) a=2;
    164164                        if (a>=0) {
    165                                 t = TriaAdjTriangles[a];
     165                                t = adj[a];
    166166                                //  if (t-this<0) return 0;
    167                                 v2 = TriaVertices[VerticesOfTriangularEdge[a][0]];
    168                                 v0 = TriaVertices[VerticesOfTriangularEdge[a][1]];
    169                                 v1 = TriaVertices[OppositeEdge[a]];
    170                                 v3 = t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
     167                                v2 = vertices[VerticesOfTriangularEdge[a][0]];
     168                                v0 = vertices[VerticesOfTriangularEdge[a][1]];
     169                                v1 = vertices[OppositeEdge[a]];
     170                                v3 = t->vertices[OppositeEdge[AdjEdgeNumber[a]&3]];
    171171                        }
    172172                }
     
    177177        double   Triangle::QualityQuad(int a,int option) const{
    178178                double q;
    179                 if (!link || TriaAdjSharedEdge[a] &4)
     179                if (!link || AdjEdgeNumber[a] &4)
    180180                 q=  -1;
    181181                else {
    182                         Triangle * t = TriaAdjTriangles[a];
     182                        Triangle * t = adj[a];
    183183                        if (t-this<0) q=  -1;// because we do 2 times
    184184                        else if (!t->link ) q=  -1;
    185                         else if (TriaAdjSharedEdge[0] & 16 || TriaAdjSharedEdge[1] & 16  || TriaAdjSharedEdge[2] & 16 || t->TriaAdjSharedEdge[0] & 16 || t->TriaAdjSharedEdge[1] & 16 || t->TriaAdjSharedEdge[2] & 16 )
     185                        else if (AdjEdgeNumber[0] & 16 || AdjEdgeNumber[1] & 16  || AdjEdgeNumber[2] & 16 || t->AdjEdgeNumber[0] & 16 || t->AdjEdgeNumber[1] & 16 || t->AdjEdgeNumber[2] & 16 )
    186186                         q= -1;
    187187                        else if(option){
    188                                 const BamgVertex & v2 = *TriaVertices[VerticesOfTriangularEdge[a][0]];
    189                                 const BamgVertex & v0 = *TriaVertices[VerticesOfTriangularEdge[a][1]];
    190                                 const BamgVertex & v1 = *TriaVertices[OppositeEdge[a]];
    191                                 const BamgVertex & v3 = * t->TriaVertices[OppositeEdge[TriaAdjSharedEdge[a]&3]];
     188                                const BamgVertex & v2 = *vertices[VerticesOfTriangularEdge[a][0]];
     189                                const BamgVertex & v0 = *vertices[VerticesOfTriangularEdge[a][1]];
     190                                const BamgVertex & v1 = *vertices[OppositeEdge[a]];
     191                                const BamgVertex & v3 = * t->vertices[OppositeEdge[AdjEdgeNumber[a]&3]];
    192192                                q =  QuadQuality(v0,v1,v2,v3); // do the float part
    193193                        }
     
    200200        void Triangle::Set(const Triangle & rec,const Mesh & Th ,Mesh & ThNew){
    201201                *this = rec;
    202                 if ( TriaVertices[0] ) TriaVertices[0] = ThNew.vertices +  Th.Number(TriaVertices[0]);
    203                 if ( TriaVertices[1] ) TriaVertices[1] = ThNew.vertices +  Th.Number(TriaVertices[1]);
    204                 if ( TriaVertices[2] ) TriaVertices[2] = ThNew.vertices +  Th.Number(TriaVertices[2]);
    205                 if(TriaAdjTriangles[0]) TriaAdjTriangles[0] =  ThNew.triangles + Th.Number(TriaAdjTriangles[0]);
    206                 if(TriaAdjTriangles[1]) TriaAdjTriangles[1] =  ThNew.triangles + Th.Number(TriaAdjTriangles[1]);
    207                 if(TriaAdjTriangles[2]) TriaAdjTriangles[2] =  ThNew.triangles + Th.Number(TriaAdjTriangles[2]);
     202                if ( vertices[0] ) vertices[0] = ThNew.vertices +  Th.Number(vertices[0]);
     203                if ( vertices[1] ) vertices[1] = ThNew.vertices +  Th.Number(vertices[1]);
     204                if ( vertices[2] ) vertices[2] = ThNew.vertices +  Th.Number(vertices[2]);
     205                if(adj[0]) adj[0] =  ThNew.triangles + Th.Number(adj[0]);
     206                if(adj[1]) adj[1] =  ThNew.triangles + Th.Number(adj[1]);
     207                if(adj[2]) adj[2] =  ThNew.triangles + Th.Number(adj[2]);
    208208                if (link  >= Th.triangles && link  < Th.triangles + Th.nbt)
    209209                 link = ThNew.triangles + Th.Number(link);
     
    216216                if(a/4 !=0) return 0;// arete lock or MarkUnSwap
    217217
    218                 register Triangle *t1=this,*t2=TriaAdjTriangles[a];// les 2 triangles adjacent
    219                 register short a1=a,a2=TriaAdjSharedEdge[a];// les 2 numero de l arete dans les 2 triangles
     218                register Triangle *t1=this,*t2=adj[a];// les 2 triangles adjacent
     219                register short a1=a,a2=AdjEdgeNumber[a];// les 2 numero de l arete dans les 2 triangles
    220220                if(a2/4 !=0) return 0; // arete lock or MarkUnSwap
    221221
    222                 register BamgVertex  *sa=t1->TriaVertices[VerticesOfTriangularEdge[a1][0]];
    223                 register BamgVertex  *sb=t1->TriaVertices[VerticesOfTriangularEdge[a1][1]];
    224                 register BamgVertex  *s1=t1->TriaVertices[OppositeVertex[a1]];
    225                 register BamgVertex  *s2=t2->TriaVertices[OppositeVertex[a2]];
     222                register BamgVertex  *sa=t1->vertices[VerticesOfTriangularEdge[a1][0]];
     223                register BamgVertex  *sb=t1->vertices[VerticesOfTriangularEdge[a1][1]];
     224                register BamgVertex  *s1=t1->vertices[OppositeVertex[a1]];
     225                register BamgVertex  *s2=t2->vertices[OppositeVertex[a2]];
    226226
    227227                Icoor2 det1=t1->det , det2=t2->det ;
  • issm/trunk/src/c/objects/Bamg/Triangle.h

    r5143 r5148  
    1717
    1818                private:
    19                         BamgVertex *TriaVertices[3];        // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
    20                         Triangle   *TriaAdjTriangles[3];    // 3 pointers toward the adjacent triangles
    21                         short       TriaAdjSharedEdge[3];   // edge id in the adjacent triangles. The edge number 1 is the edge number TriaAdjSharedEdge[1] in the Adjacent triangle 1
     19                        BamgVertex *vertices[3];        // 3 vertices if t is triangle, t[i] allowed by access function, (*t)[i] if pointer
     20                        Triangle   *adj[3];    // 3 pointers toward the adjacent triangles
     21                        short       AdjEdgeNumber[3];   // edge id in the adjacent triangles. The edge number 1 is the edge number AdjEdgeNumber[1] in the Adjacent triangle 1
    2222
    2323                public:
     
    3434
    3535                        //Operators
    36                         const BamgVertex & operator[](int i) const {return *TriaVertices[i];};
    37                         BamgVertex & operator[](int i)  {return *TriaVertices[i];};
    38                         const BamgVertex * operator()(int i) const {return TriaVertices[i];};
    39                         BamgVertex * & operator()(int i)  {return TriaVertices[i];};
     36                        const BamgVertex & operator[](int i) const {return *vertices[i];};
     37                        BamgVertex & operator[](int i)  {return *vertices[i];};
     38                        const BamgVertex * operator()(int i) const {return vertices[i];};
     39                        BamgVertex * & operator()(int i)  {return vertices[i];};
    4040
    4141                        //Methods
     
    4343                        int    swap(short a1,int=0);
    4444                        long   Optim(short a,int =0);
    45                         int    Locked(int a)const { return TriaAdjSharedEdge[a]&4;}
    46                         int    Hidden(int a)const { return TriaAdjSharedEdge[a]&16;}
    47                         int    GetAllflag(int a){return TriaAdjSharedEdge[a] & 1020;}
    48                         void   SetAllFlag(int a,int f){TriaAdjSharedEdge[a] = (TriaAdjSharedEdge[a] &3) + (1020 & f);}
     45                        int    Locked(int a)const { return AdjEdgeNumber[a]&4;}
     46                        int    Hidden(int a)const { return AdjEdgeNumber[a]&16;}
     47                        int    GetAllflag(int a){return AdjEdgeNumber[a] & 1020;}
     48                        void   SetAllFlag(int a,int f){AdjEdgeNumber[a] = (AdjEdgeNumber[a] &3) + (1020 & f);}
    4949                        double QualityQuad(int a,int option=1) const;
    50                         short  NuEdgeTriangleAdj(int i) const {return TriaAdjSharedEdge[i&3]&3;} // Number of the  adjacent edge in adj tria 
     50                        short  NuEdgeTriangleAdj(int i) const {return AdjEdgeNumber[i&3]&3;} // Number of the  adjacent edge in adj tria 
    5151                        AdjacentTriangle FindBoundaryEdge(int i) const;
    52                         AdjacentTriangle Adj(int i)  const {return AdjacentTriangle(TriaAdjTriangles[i],TriaAdjSharedEdge[i]&3);};
    53                         Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}
     52                        AdjacentTriangle Adj(int i)  const {return AdjacentTriangle(adj[i],AdjEdgeNumber[i]&3);};
     53                        Triangle* TriangleAdj(int i) const {return adj[i&3];}
    5454                        Triangle* Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const ;
    5555                        void  Renumbering(Triangle *tb,Triangle *te, long *renu){
    5656                                if (link  >=tb && link  <te) link  = tb + renu[link -tb];
    57                                 if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb];
    58                                 if (TriaAdjTriangles[1] >=tb && TriaAdjTriangles[1] <te) TriaAdjTriangles[1] = tb + renu[TriaAdjTriangles[1]-tb];
    59                                 if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];   
     57                                if (adj[0] >=tb && adj[0] <te) adj[0] = tb + renu[adj[0]-tb];
     58                                if (adj[1] >=tb && adj[1] <te) adj[1] = tb + renu[adj[1]-tb];
     59                                if (adj[2] >=tb && adj[2] <te) adj[2] = tb + renu[adj[2]-tb];   
    6060                        }
    6161                        void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
    62                                 if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
    63                                 if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
    64                                 if (TriaVertices[2] >=vb && TriaVertices[2] <ve) TriaVertices[2] = vb + renu[TriaVertices[2]-vb];   
     62                                if (vertices[0] >=vb && vertices[0] <ve) vertices[0] = vb + renu[vertices[0]-vb];
     63                                if (vertices[1] >=vb && vertices[1] <ve) vertices[1] = vb + renu[vertices[1]-vb];
     64                                if (vertices[2] >=vb && vertices[2] <ve) vertices[2] = vb + renu[vertices[2]-vb];   
    6565                        }
    6666                        void SetAdjAdj(short a){
    6767                                a &= 3;
    68                                 register Triangle *tt=TriaAdjTriangles[a];
    69                                 TriaAdjSharedEdge [a] &= 55; // remove MarkUnSwap
    70                                 register short aatt = TriaAdjSharedEdge[a] & 3;
     68                                register Triangle *tt=adj[a];
     69                                AdjEdgeNumber [a] &= 55; // remove MarkUnSwap
     70                                register short aatt = AdjEdgeNumber[a] & 3;
    7171                                if(tt){
    72                                         tt->TriaAdjTriangles[aatt]=this;
    73                                         tt->TriaAdjSharedEdge[aatt]=a + (TriaAdjSharedEdge[a] & 60 ) ;}// Copy all the mark
     72                                        tt->adj[aatt]=this;
     73                                        tt->AdjEdgeNumber[aatt]=a + (AdjEdgeNumber[a] & 60 ) ;}// Copy all the mark
    7474                          }
    7575                        void SetAdj2(short a,Triangle *t,short aat){
    76                                 TriaAdjTriangles[a]=t;    //the adjacent triangle to the edge a is t
    77                                 TriaAdjSharedEdge[a]=aat; //position of the edge in the adjacent triangle
     76                                adj[a]=t;    //the adjacent triangle to the edge a is t
     77                                AdjEdgeNumber[a]=aat; //position of the edge in the adjacent triangle
    7878                                if(t) { //if t!=NULL add adjacent triangle to t (this)
    79                                         t->TriaAdjTriangles[aat]=this;
    80                                         t->TriaAdjSharedEdge[aat]=a;
     79                                        t->adj[aat]=this;
     80                                        t->AdjEdgeNumber[aat]=a;
    8181                                }
    8282                        }
    8383                        void SetSingleVertexToTriangleConnectivity() {
    84                                 if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0);
    85                                 if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1);
    86                                 if (TriaVertices[2]) (TriaVertices[2]->t=this,TriaVertices[2]->vint=2);
     84                                if (vertices[0]) (vertices[0]->t=this,vertices[0]->vint=0);
     85                                if (vertices[1]) (vertices[1]->t=this,vertices[1]->vint=1);
     86                                if (vertices[2]) (vertices[2]->t=this,vertices[2]->vint=2);
    8787                        }
    8888                        void SetHidden(int a){
    8989                                //Get Adjacent Triangle number a
    90                                 register Triangle* t = TriaAdjTriangles[a];
     90                                register Triangle* t = adj[a];
    9191                                //if it exist
    9292                                //C|=D -> C=(C|D) bitwise inclusive OR
    93                                 if(t) t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=16;
    94                                 TriaAdjSharedEdge[a] |= 16;
     93                                if(t) t->AdjEdgeNumber[AdjEdgeNumber[a] & 3] |=16;
     94                                AdjEdgeNumber[a] |= 16;
    9595                        }
    9696
    9797                        void SetLocked(int a){
    9898                                //mark the edge as on Boundary
    99                                 register Triangle * t = TriaAdjTriangles[a];
    100                                 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=4;
    101                                 TriaAdjSharedEdge[a] |= 4;
     99                                register Triangle * t = adj[a];
     100                                t->AdjEdgeNumber[AdjEdgeNumber[a] & 3] |=4;
     101                                AdjEdgeNumber[a] |= 4;
    102102                        }
    103103                        void SetMarkUnSwap(int a){
    104                                 register Triangle * t = TriaAdjTriangles[a];
    105                                 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] |=8;
    106                                 TriaAdjSharedEdge[a] |=8 ;
     104                                register Triangle * t = adj[a];
     105                                t->AdjEdgeNumber[AdjEdgeNumber[a] & 3] |=8;
     106                                AdjEdgeNumber[a] |=8 ;
    107107                        }
    108108                        void SetUnMarkUnSwap(int a){
    109                                 register Triangle * t = TriaAdjTriangles[a];
    110                                 t->TriaAdjSharedEdge[TriaAdjSharedEdge[a] & 3] &=55; // 23 + 32
    111                                 TriaAdjSharedEdge[a] &=55 ;
     109                                register Triangle * t = adj[a];
     110                                t->AdjEdgeNumber[AdjEdgeNumber[a] & 3] &=55; // 23 + 32
     111                                AdjEdgeNumber[a] &=55 ;
    112112                        }
    113113                        void SetDet() {
    114                                 if(TriaVertices[0] && TriaVertices[1] && TriaVertices[2])    det = bamg::det(*TriaVertices[0],*TriaVertices[1],*TriaVertices[2]);
     114                                if(vertices[0] && vertices[1] && vertices[2])    det = bamg::det(*vertices[0],*vertices[1],*vertices[2]);
    115115                                else det = -1; }
    116116
     
    118118                        double qualite() ;
    119119                        void  Set(const Triangle &,const Mesh &,Mesh &);
    120                         int   In(BamgVertex *v) const { return TriaVertices[0]==v || TriaVertices[1]==v || TriaVertices[2]==v ;}
     120                        int   In(BamgVertex *v) const { return vertices[0]==v || vertices[1]==v || vertices[2]==v ;}
    121121
    122122        };
  • issm/trunk/src/c/objects/Bamg/VertexOnEdge.h

    r5120 r5148  
    2929
    3030                        //Methods
    31                         void SetOnBTh(){v->onBackgroundEdge=this;v->vint=IsVertexOnEdge;} 
     31                        void SetOnBTh(){v->BackgroundEdgeHook=this;v->vint=IsVertexOnEdge;} 
    3232                        void Set(const Mesh &,long,Mesh &); 
    3333        };
  • issm/trunk/src/c/objects/Bamg/VertexOnGeom.h

    r5120 r5148  
    3838                        int  OnGeomEdge() const {return this? abscisse >=0 :0;}
    3939                        int  IsRequiredVertex() {return this? ((abscisse<0 ? (gv?gv->Required():0):0 )) : 0;}
    40                         void SetOn(){mv->onGeometry=this;mv->vint=IsVertexOnGeom;}
     40                        void SetOn(){mv->GeometricalEdgeHook=this;mv->vint=IsVertexOnGeom;}
    4141
    4242                        //Inline methods
  • issm/trunk/src/c/objects/Bamg/VertexOnVertex.h

    r5120 r5148  
    2121
    2222                        //Methods
    23                         void SetOnBTh(){v->onBackgroundVertex=bv;v->vint=IsVertexOnVertex;}
     23                        void SetOnBTh(){v->BackgroundVertexHook=bv;v->vint=IsVertexOnVertex;}
    2424                        void Set(const Mesh &,long,Mesh &);
    2525        };
Note: See TracChangeset for help on using the changeset viewer.