Changeset 5124


Ignore:
Timestamp:
08/10/10 14:26:35 (15 years ago)
Author:
Mathieu Morlighem
Message:

Renamed several Bamg routines

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

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/modules/Bamgx/Bamgx.cpp

    r5095 r5124  
    6767
    6868                //Renumbering
    69                 Th.ReNumberingTheTriangleBySubDomain();
     69                Th.TrianglesRenumberBySubDomain();
    7070
    7171                //Crack mesh if requested
    … …  
    177177
    178178                //Renumber by subdomain
    179                 Th.ReNumberingTheTriangleBySubDomain();
     179                Th.TrianglesRenumberBySubDomain();
    180180
    181181                //Smooth vertices

  • issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp

    r5095 r5124  
    6161        if (verbose) printf("Reading mesh\n");
    6262        Mesh Th(index_data,x_data,y_data,nods_data,nels_data);
    63         Th.ReMakeTriangleContainingTheVertex();
     63        Th.CreateSingleVertexToTriangleConnectivity();
    6464
    6565        //Loop over output nodes
    … …  
    8585
    8686                        //Find triangle holding r/I
    87                         Triangle &tb=*Th.FindTriangleContaining(I,dete);
     87                        Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
    8888
    8989                        // internal point

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

    r5120 r5124  
    4848                I2 IBTh  = BTh.toI2(PNew);
    4949
    50                 tstart=BTh.FindTriangleContaining(IBTh,deta,tstart); 
     50                tstart=BTh.TriangleFindFromCoord(IBTh,deta,tstart); 
    5151
    5252                if (tstart->det <0){ // outside

  • issm/trunk/src/c/objects/Bamg/Edge.h

    r5120 r5124  
    2828
    2929                        //Methods
    30                         void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
     30                        void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
    3131                                if (v[0] >=vb && v[0] <ve) v[0] = vb + renu[v[0]-vb];
    3232                                if (v[1] >=vb && v[1] <ve) v[1] = vb + renu[v[1]-vb];

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

    r3913 r5124  
    1616        /*FUNCTION GeometricalVertex::Corner {{{1*/
    1717        int  GeometricalVertex::Corner() const {
    18                 return cas & 4;
     18                return type & 4;
    1919        }
    2020        /*}}}*/
    … …  
    2222        int  GeometricalVertex::Required()const {
    2323                // a corner is required
    24                 return cas & 6;
     24                return type & 6;
    2525        }
    26         /*}}}*/
    27         /*FUNCTION GeometricalVertex::IsThe {{{1*/
    28         int  GeometricalVertex::IsThe() const {
    29                 return link == this;
    30         } 
    3126        /*}}}*/
    3227        /*FUNCTION GeometricalVertex::SetCorner {{{1*/
    3328        void GeometricalVertex::SetCorner(){
    34                 cas |= 4;
     29                type |= 4;
    3530        }
    3631        /*}}}*/
    3732        /*FUNCTION GeometricalVertex::SetRequired {{{1*/
    3833        void GeometricalVertex::SetRequired(){
    39                 cas |= 2;
     34                type |= 2;
    4035        }
    4136        /*}}}*/
    42         /*FUNCTION GeometricalVertex::Set {{{1*/
    43         void GeometricalVertex::Set(){
    44                 cas=0;
    45         }
    46         /*}}}*/
    47         /*FUNCTION GeometricalVertex::The {{{1*/
    48         GeometricalVertex* GeometricalVertex::The(){
    49                 // return a unique vertex
    50                 ISSMASSERT(link);
    51                 return link;
    52         }/*}}}*/
    5337
    5438}

  • issm/trunk/src/c/objects/Bamg/GeometricalVertex.h

    r5120 r5124  
    1414                        friend class Geometry;
    1515
    16                         int cas;
    17                         GeometricalVertex* link; //  link all the same GeometricalVertex circular (Crack)
     16                        int type;
    1817
    1918                        //Constructors
    20                         GeometricalVertex() :cas(0), link(this) {};
     19                        GeometricalVertex():type(0){};
    2120
    2221                        //Methods
    2322                        int  Corner() const;
    2423                        int  Required()const;
    25                         int  IsThe() const;
    2624                        void SetCorner();
    2725                        void SetRequired();
    28                         void Set();
    29                         GeometricalVertex* The();
    3026
    31                         //Inline methods
    32                         inline void Set(const GeometricalVertex & rec,const Geometry & ,const Geometry & ){
    33                                 *this=rec;
    34                         }
    3527        };
    3628

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

    r5120 r5124  
    3333                curves= NbOfCurves ? new Curve[NbOfCurves]:NULL;
    3434                subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL;
    35                 for (i=0;i<nbv;i++)
    36                  vertices[i].Set(Gh.vertices[i],Gh,*this);
    3735                for (i=0;i<nbe;i++)
    3836                 edges[i].Set(Gh.edges[i],Gh,*this);
    … …  
    8987                                vertices[i].DirOfSearch=NoDirOfSearch;
    9088                                vertices[i].color =0;
    91                                 vertices[i].Set();
     89                                vertices[i].type=0;
    9290                        }
    9391                        /*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/
    … …  
    424422
    425423                k=0;
    426                 //link all vertices to themselves by default
    427                 for (i=0;i<nbv;i++) vertices[i].link = vertices+i;
    428424
    429425                //build quadtree for this geometry
    … …  
    443439                        //if there is a vertex found that is to close to vertices[i] -> error
    444440                        if( v && Norme1(v->r - vertices[i]) < eps ){
    445                                 // mama's old trick to get j
    446                                 GeometricalVertex* vg=(GeometricalVertex*) (void*)v;
    447                                 j=vg-v0g;
    448                                 //check that the clostest vertex is not itself...
    449                                 if ( v !=  &(BamgVertex &) vertices[j]){
    450                                         ISSMERROR(" v !=  &(BamgVertex &) vertices[j]");
    451                                 }
    452                                 vertices[i].link = vertices + j;
    453                                 k++;         
    454                         }
    455 
    456                         //The nearest vertex was non existent or far enough from vertices[i]
     441                                printf("WARNING: two points of the geometry are very closed to each other\n");
     442                        }
     443
    457444                        //Add vertices[i] to the quadtree
    458                         else  quadtree.Add(vertices[i]);
     445                        quadtree.Add(vertices[i]);
    459446                }
    460447
    … …  
    463450                        printf("number of distinct vertices= %i, over %i\n",nbv - k,nbv);
    464451                        printf("List of duplicate vertices:\n");
    465                         for (i=0;i<nbv;i++){
    466                                 if (!vertices[i].IsThe()) printf("  %i and %i\n",i,Number(vertices[i].The()));
    467                         }
    468452                        ISSMERROR("See above");
    469453                }

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

    r5120 r5124  
    3434                else {// not optimisation
    3535                        init();
    36                         t=tbegin = Bh.FindTriangleContaining(a,deta);
     36                        t=tbegin = Bh.TriangleFindFromCoord(a,deta);
    3737                        if( t->det>=0)
    3838                         ilast=NewItem(t,double(deta[0])/t->det,double(deta[1])/t->det,double(deta[2])/t->det);

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

    r5120 r5124  
    11881188                tt[i2]->SetAdj2(i1,tt[i1],i2);
    11891189
    1190                 tt[0]->SetTriangleContainingTheVertex();
    1191                 tt[1]->SetTriangleContainingTheVertex();
    1192                 tt[2]->SetTriangleContainingTheVertex();
     1190                tt[0]->SetSingleVertexToTriangleConnectivity();
     1191                tt[1]->SetSingleVertexToTriangleConnectivity();
     1192                tt[2]->SetSingleVertexToTriangleConnectivity();
    11931193
    11941194
    … …  
    22562256
    22572257                //Now, find the triangles that hold a cracked edge
    2258                 ReMakeTriangleContainingTheVertex();
     2258                CreateSingleVertexToTriangleConnectivity();
    22592259
    22602260                long* Edgeflags=new long[NbCrackedEdges];
    … …  
    25372537                                NbSubDomains =Gh.NbSubDomains;
    25382538                                long err=0;
    2539                                 ReMakeTriangleContainingTheVertex();
     2539                                CreateSingleVertexToTriangleConnectivity();
    25402540                                long * mark = new long[nbt];
    25412541                                Edge **GeometricalEdgetoEdge = MakeGeometricalEdgeToEdge();
    … …  
    26232623        }
    26242624        /*}}}1*/
    2625         /*FUNCTION Mesh::FindTriangleContaining{{{1*/
    2626         Triangle * Mesh::FindTriangleContaining(const I2 & B,Icoor2 dete[3], Triangle *tstart) const {
    2627                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/
    2628 
    2629                 Triangle * t=0;
    2630                 int j,jp,jn,jj;
    2631                 int counter;
    2632 
    2633                 /*Get starting triangle. Take tsart if provided*/
    2634                 if (tstart){
    2635                         t=tstart;
    2636                 }
    2637                 /*Else find the closest Triangle using the quadtree*/
    2638                 else {
    2639 
    2640                         /*Check that the quadtree does exist*/
    2641                         if (!quadtree) ISSMERROR("no starting triangle provided and no quadtree available");
    2642 
    2643                         /*Call NearestVertex*/
    2644                         BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
    2645 
    2646                         /*Check output (Vertex a)*/
    2647                         if (!a)    ISSMERROR("problem while trying to find nearest vertex from a given point. No output found");
    2648                         if (!a->t) ISSMERROR("no triangle is associated to vertex number %i (another call to ReMakeTriangleContainingTheVertex is required)",Number(a)+1);
    2649                         ISSMASSERT(a>=vertices && a<vertices+nbv);
    2650 
    2651                         /*Get starting triangle*/
    2652                         t = a->t;
    2653                         ISSMASSERT(t>=triangles && t<triangles+nbt);
    2654                 }
    2655 
    2656                 Icoor2  detop ;
    2657 
    2658                 /*initialize number of test triangle*/
    2659                 counter=0;
    2660 
    2661                 /*The initial triangle might be outside*/
    2662                 while (t->det < 0){
    2663 
    2664                         /*Get a real vertex from this triangle (k0)*/
    2665                         int k0=(*t)(0)?(((*t)(1)?((*t)(2)?-1:2):1)):0;
    2666                         ISSMASSERT(k0>=0);// k0 the NULL vertex
    2667                         int k1=NextVertex[k0],k2=PreviousVertex[k0];
    2668                         dete[k0]=det(B,(*t)[k1],(*t)[k2]);
    2669                         dete[k1]=dete[k2]=-1;     
    2670                         if (dete[k0] > 0) // outside B
    2671                          return t;
    2672                         t = t->TriangleAdj(OppositeEdge[k0]);
    2673                         counter++;
    2674                         ISSMASSERT(counter<2);
    2675                 }
    2676 
    2677                 jj=0;
    2678                 detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
    2679 
    2680                 while(t->det>0) {
    2681 
    2682                         /*Increase counter*/
    2683                         if (++counter>=10000) ISSMERROR("Maximum number of iteration reached (threshold = %i).",counter);
    2684 
    2685                         j= OppositeVertex[jj];
    2686                         dete[j] = detop;  //det(*b,*s1,*s2);
    2687                         jn = NextVertex[j];
    2688                         jp = PreviousVertex[j];
    2689                         dete[jp]= det(*(*t)(j),*(*t)(jn),B);
    2690                         dete[jn] = t->det-dete[j] -dete[jp];
    2691 
    2692                         // count the number k of  dete <0
    2693                         int k=0,ii[3];
    2694                         if (dete[0] < 0 ) ii[k++]=0;
    2695                         if (dete[1] < 0 ) ii[k++]=1;
    2696                         if (dete[2] < 0 ) ii[k++]=2;
    2697                         // 0 => ok
    2698                         // 1 => go in way 1
    2699                         // 2 => two way go in way 1 or 2 randomly
    2700 
    2701                         if (k==0) break;
    2702                         if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]);
    2703                         if ( k>=3){
    2704                                 ISSMERROR("k>=3");
    2705                         }
    2706                         TriangleAdjacent t1 = t->Adj(jj=ii[0]);
    2707                         if ((t1.det() < 0 ) && (k == 2))
    2708                          t1 = t->Adj(jj=ii[1]);
    2709                         t=t1;
    2710                         j=t1;// for optimisation we now the -det[OppositeVertex[j]];
    2711                         detop = -dete[OppositeVertex[jj]];
    2712                         jj = j;
    2713                 }
    2714 
    2715                 if (t->det<0) // outside triangle
    2716                  dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;
    2717                 //  NbOfTriangleSearchFind += counter; 
    2718                 return t;
    2719         }
    2720         /*}}}1*/
    27212625        /*FUNCTION Mesh::Init{{{1*/
    27222626        void Mesh::Init(long maxnbv_in) {
    … …  
    28612765                triangles[1].det = -1;  //boundary triangle: det = -1
    28622766
    2863                 triangles[0].SetTriangleContainingTheVertex();
    2864                 triangles[1].SetTriangleContainingTheVertex();
     2767                triangles[0].SetSingleVertexToTriangleConnectivity();
     2768                triangles[1].SetSingleVertexToTriangleConnectivity();
    28652769
    28662770                triangles[0].link=&triangles[1];
    … …  
    28842788                        //Find the triangle in which newvertex is located
    28852789                        Icoor2 dete[3];
    2886                         Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates)
     2790                        Triangle* tcvi = TriangleFindFromCoord(newvertex->i,dete); //(newvertex->i = integer coordinates)
    28872791
    28882792                        //Add newvertex to the quadtree
    … …  
    29182822                        if(!NbSwap) break;
    29192823                }
    2920                 ReMakeTriangleContainingTheVertex();
     2824                CreateSingleVertexToTriangleConnectivity();
    29212825                // because we break the TriangleContainingTheVertex
    29222826#endif
    … …  
    29752879                                }
    29762880                                vj.ReferenceNumber=0;
    2977                                 Triangle *tcvj=FindTriangleContaining(vj.i,dete);
     2881                                Triangle *tcvj=TriangleFindFromCoord(vj.i,dete);
    29782882                                if (tcvj && !tcvj->link){
    29792883                                        tcvj->Echo();
    … …  
    31553059                I2 a = toI2(A);
    31563060                Icoor2 deta[3];
    3157                 Triangle * t =FindTriangleContaining(a,deta);
     3061                Triangle * t =TriangleFindFromCoord(a,deta);
    31583062                if (t->det <0) { // outside
    31593063                        double ba,bb;
    … …  
    32023106                                }
    32033107                        }
    3204                         Bh.ReMakeTriangleContainingTheVertex();     
     3108                        Bh.CreateSingleVertexToTriangleConnectivity();     
    32053109                        InsertNewPoints(nbvold,NbTSwap)   ;           
    32063110                } 
    3207                 else Bh.ReMakeTriangleContainingTheVertex();     
     3111                else Bh.CreateSingleVertexToTriangleConnectivity();     
    32083112
    32093113                // generation of the list of next Triangle
    … …  
    33703274                double abscisse = -1;
    33713275
    3372                 for (int cas=0;cas<2;cas++){
     3276                for (int step=0;step<2;step++){
    33733277                        // 2 times algo:
    33743278                        //    1 for computing the length l
    … …  
    33853289
    33863290                                kkk=kkk+1;
    3387                                 if (kkk>=100){
    3388                                         ISSMERROR("kkk>=100");
    3389                                 }
    3390                                 if (!eee){
    3391                                         ISSMERROR("!eee");
    3392                                 }
     3291                                ISSMASSERT(kkk<100);
     3292                                ISSMASSERT(eee);
    33933293                                double lg0 = lg;
    33943294                                double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
    33953295                                lg += dp;
    3396                                 if (cas && abscisse <= lg) { // ok we find the geom edge
     3296                                if (step && abscisse <= lg) { // ok we find the geom edge
    33973297                                        double sss  =   (abscisse-lg0)/dp;
    33983298                                        double thetab = te0*(1-sss)+ sss*iii;
    3399                                         if (thetab<0 || thetab>1){
    3400                                                 ISSMERROR("thetab<0 || thetab>1");
    3401                                         }
     3299                                        ISSMASSERT(thetab>=0 && thetab<=1);
    34023300                                        BR = VertexOnEdge(&R,eee,thetab);
    34033301                                        return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
    … …  
    34103308
    34113309                                double lg0 = lg;
    3412                                 if (!eee){
    3413                                         ISSMERROR("!eee");
    3414                                 }
     3310                                ISSMASSERT(eee);
    34153311                                v1 = pvB;
    34163312                                double dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0);
    … …  
    34213317                                        double sss  =   (abscisse-lg0)/dp;
    34223318                                        double thetab = te0*(1-sss)+ sss*tB;
    3423                                         if (thetab<0 || thetab>1){
    3424                                                 ISSMERROR("thetab<0 || thetab>1");
    3425                                         }
     3319                                        ISSMASSERT(thetab>=0 && thetab<=1);
    34263320                                        BR = VertexOnEdge(&R,eee,thetab);
    34273321                                        return  Gh.ProjectOnCurve(*eee,thetab,R,GR);
    … …  
    35943488        triangles[1].det = -1;  // boundary triangles
    35953489
    3596         triangles[0].SetTriangleContainingTheVertex();
    3597         triangles[1].SetTriangleContainingTheVertex();
     3490        triangles[0].SetSingleVertexToTriangleConnectivity();
     3491        triangles[1].SetSingleVertexToTriangleConnectivity();
    35983492
    35993493        triangles[0].link=&triangles[1];
    … …  
    36103504                BamgVertex *vi  = ordre[icount];
    36113505                Icoor2 dete[3];
    3612                 Triangle *tcvi = FindTriangleContaining(vi->i,dete);
     3506                Triangle *tcvi = TriangleFindFromCoord(vi->i,dete);
    36133507                quadtree->Add(*vi);
    36143508                AddVertex(*vi,tcvi,dete);
    … …  
    37423636}
    37433637/*}}}1*/
    3744         /*FUNCTION Mesh::ReNumberingTheTriangleBySubDomain{{{1*/
    3745         void Mesh::ReNumberingTheTriangleBySubDomain(bool justcompress){
     3638        /*FUNCTION Mesh::TrianglesRenumberBySubDomain{{{1*/
     3639        void Mesh::TrianglesRenumberBySubDomain(bool justcompress){
    37463640                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/
    37473641
    … …  
    37863680                // do the change on all the pointeur
    37873681                for ( it=0;it<nbt;it++)
    3788                  triangles[it].ReNumbering(triangles,te,renu);
     3682                 triangles[it].Renumbering(triangles,te,renu);
    37893683
    37903684                for ( i=0;i<NbSubDomains;i++)
    … …  
    38123706        }
    38133707        /*}}}1*/
    3814         /*FUNCTION Mesh::ReNumberingVertex{{{1*/
    3815         void Mesh::ReNumberingVertex(long * renu) {
     3708        /*FUNCTION Mesh::VerticesRenumber{{{1*/
     3709        void Mesh::VerticesRenumber(long * renu) {
    38163710                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingVertex)*/
    38173711
    … …  
    38243718                printf("renumbering triangles\n");
    38253719                for ( it=0;it<nbt;it++)
    3826                  triangles[it].ReNumbering(vertices,ve,renu);
     3720                 triangles[it].Renumbering(vertices,ve,renu);
    38273721
    38283722                printf("renumbering edges\n");
    38293723                for ( ie=0;ie<nbe;ie++)
    3830                  edges[ie].ReNumbering(vertices,ve,renu);
     3724                 edges[ie].Renumbering(vertices,ve,renu);
    38313725
    38323726                printf("renumbering vertices on geom\n");
    … …  
    40673961        if (quadtree) delete quadtree;
    40683962        quadtree=0;
    4069         ReMakeTriangleContainingTheVertex();
     3963        CreateSingleVertexToTriangleConnectivity();
    40703964        Triangle vide; // a triangle to mark the boundary vertex
    40713965        Triangle   ** tstart= new Triangle* [nbv];
    … …  
    41094003        if(raisonmax<1.1) return;
    41104004        if(verbose > 1) printf("   Mesh::SmoothMetric raisonmax = %g\n",raisonmax);
    4111         ReMakeTriangleContainingTheVertex();
     4005        CreateSingleVertexToTriangleConnectivity();
    41124006        long i,j,kch,kk,ip;
    41134007        long *first_np_or_next_t0 = new long[nbv];
    … …  
    41884082                const  int withBackground = &BTh != this && &BTh;
    41894083
    4190                 ReNumberingTheTriangleBySubDomain();
     4084                TrianglesRenumberBySubDomain();
    41914085                //int nswap =0;
    41924086                const long nfortria( choice ? 4 : 6);
    … …  
    46384532                SetIntCoor("In SplitElement");
    46394533
    4640                 ReMakeTriangleContainingTheVertex();
     4534                CreateSingleVertexToTriangleConnectivity();
    46414535                if(withBackground)
    4642                  BTh.ReMakeTriangleContainingTheVertex();
     4536                 BTh.CreateSingleVertexToTriangleConnectivity();
    46434537
    46444538                delete [] edges;
    … …  
    46624556                NbVerticesOnGeomEdge = newNbVerticesOnGeomEdge;
    46634557                NbVertexOnBThEdge=newNbVertexOnBThEdge;
    4664                 //  ReMakeTriangleContainingTheVertex();
     4558                //  CreateSingleVertexToTriangleConnectivity();
    46654559
    46664560                ReconstructExistingMesh();
    … …  
    47214615                  }
    47224616        }
    4723         ReMakeTriangleContainingTheVertex();   
     4617        CreateSingleVertexToTriangleConnectivity();   
    47244618        if (nbvold!=nbv){
    47254619                long  iv = nbvold;
    … …  
    47344628                        vi.ReferenceNumber=0;
    47354629                        vi.DirOfSearch =NoDirOfSearch;
    4736                         Triangle *tcvi = FindTriangleContaining(vi.i,dete);
     4630                        Triangle *tcvi = TriangleFindFromCoord(vi.i,dete);
    47374631                        if (tcvi && !tcvi->link) {
    47384632                                printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
    … …  
    47584652
    47594653return  NbSplitEdge;
     4654}
     4655/*}}}1*/
     4656/*FUNCTION Mesh::TriangleFindFromCoord{{{1*/
     4657Triangle * Mesh::TriangleFindFromCoord(const I2 & B,Icoor2 dete[3], Triangle *tstart) const {
     4658        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/
     4659
     4660        Triangle * t=0;
     4661        int j,jp,jn,jj;
     4662        int counter;
     4663
     4664        /*Get starting triangle. Take tsart if provided*/
     4665        if (tstart) t=tstart;
     4666
     4667        /*Else find the closest Triangle using the quadtree*/
     4668        else {
     4669
     4670                /*Check that the quadtree does exist*/
     4671                if (!quadtree) ISSMERROR("no starting triangle provided and no quadtree available");
     4672
     4673                /*Call NearestVertex*/
     4674                BamgVertex *a = quadtree->NearestVertex(B.x,B.y) ;
     4675
     4676                /*Check output (Vertex a)*/
     4677                if (!a)    ISSMERROR("problem while trying to find nearest vertex from a given point. No output found");
     4678                if (!a->t) ISSMERROR("no triangle is associated to vertex number %i (orphan?)",Number(a)+1);
     4679                ISSMASSERT(a>=vertices && a<vertices+nbv);
     4680
     4681                /*Get starting triangle*/
     4682                t = a->t;
     4683                ISSMASSERT(t>=triangles && t<triangles+nbt);
     4684        }
     4685
     4686        Icoor2  detop ;
     4687
     4688        /*initialize number of test triangle*/
     4689        counter=0;
     4690
     4691        /*The initial triangle might be outside*/
     4692        while (t->det < 0){
     4693
     4694                /*Get a real vertex from this triangle (k0)*/
     4695                int k0=(*t)(0)?(((*t)(1)?((*t)(2)?-1:2):1)):0;
     4696                ISSMASSERT(k0>=0);// k0 the NULL vertex
     4697                int k1=NextVertex[k0],k2=PreviousVertex[k0];
     4698                dete[k0]=det(B,(*t)[k1],(*t)[k2]);
     4699                dete[k1]=dete[k2]=-1;     
     4700                if (dete[k0] > 0) // outside B
     4701                 return t;
     4702                t = t->TriangleAdj(OppositeEdge[k0]);
     4703                counter++;
     4704                ISSMASSERT(counter<2);
     4705        }
     4706
     4707        jj=0;
     4708        detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B);
     4709
     4710        while(t->det>0) {
     4711
     4712                /*Increase counter*/
     4713                if (++counter>=10000) ISSMERROR("Maximum number of iteration reached (threshold = %i).",counter);
     4714
     4715                j= OppositeVertex[jj];
     4716                dete[j] = detop;  //det(*b,*s1,*s2);
     4717                jn = NextVertex[j];
     4718                jp = PreviousVertex[j];
     4719                dete[jp]= det(*(*t)(j),*(*t)(jn),B);
     4720                dete[jn] = t->det-dete[j] -dete[jp];
     4721
     4722                // count the number k of  dete <0
     4723                int k=0,ii[3];
     4724                if (dete[0] < 0 ) ii[k++]=0;
     4725                if (dete[1] < 0 ) ii[k++]=1;
     4726                if (dete[2] < 0 ) ii[k++]=2;
     4727                // 0 => ok
     4728                // 1 => go in way 1
     4729                // 2 => two way go in way 1 or 2 randomly
     4730
     4731                if (k==0) break;
     4732                if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]);
     4733                if ( k>=3){
     4734                        ISSMERROR("k>=3");
     4735                }
     4736                TriangleAdjacent t1 = t->Adj(jj=ii[0]);
     4737                if ((t1.det() < 0 ) && (k == 2))
     4738                 t1 = t->Adj(jj=ii[1]);
     4739                t=t1;
     4740                j=t1;// for optimisation we now the -det[OppositeVertex[j]];
     4741                detop = -dete[OppositeVertex[jj]];
     4742                jj = j;
     4743        }
     4744
     4745        if (t->det<0) // outside triangle
     4746         dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop;
     4747        //  NbOfTriangleSearchFind += counter; 
     4748        return t;
     4749}
     4750/*}}}1*/
     4751/*FUNCTION Mesh::TriangleIntNumbering{{{1*/
     4752void Mesh::TriangleIntNumbering(long* renumbering){
     4753
     4754        long num=0;
     4755        for (int i=0;i<nbt;i++){
     4756                if (triangles[i].det>0) renumbering[i]=num++;
     4757                else renumbering[i]=-1;
     4758        }
     4759        return;   
     4760}
     4761/*}}}1*/
     4762/*FUNCTION Mesh::TriangleReferenceList{{{1*/
     4763long  Mesh::TriangleReferenceList(long* reft) const {
     4764        /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)*/
     4765
     4766        long int verbose=0;
     4767        register Triangle *t0,*t;
     4768        register long k=0, num;   
     4769
     4770        //initialize all triangles as -1 (outside)
     4771        for (int it=0;it<nbt;it++) reft[it]=-1;
     4772
     4773        //loop over all subdomains
     4774        for (int i=0;i<NbSubDomains;i++){
     4775
     4776                //first triangle of the subdomain i
     4777                t=t0=subdomains[i].head;
     4778
     4779                //check that the subdomain is not empty
     4780                if (!t0){ ISSMERROR("At least one subdomain is empty");}
     4781
     4782                //loop
     4783                do{
     4784                        k++;
     4785
     4786                        //get current triangle number
     4787                        num = Number(t);
     4788
     4789                        //check that num is in [0 nbt[
     4790                        if (num<0 || num>=nbt){ ISSMERROR("num<0 || num>=nbt");}
     4791
     4792                        //reft of this triangle is the subdomain number
     4793                        reft[num]=i;
     4794
     4795                } while (t0 != (t=t->link));
     4796                //stop when all triangles of subdomains have been tagged
     4797
     4798        }
     4799        return k;   
    47604800}
    47614801/*}}}1*/
    … …  
    47934833                //Compute number of vertices on geometrical vertex
    47944834                for (i=0;i<Gh.nbv;i++){
    4795                         if (Gh[i].Required() && Gh[i].IsThe()) NbVerticesOnGeomVertex++;
     4835                        if (Gh[i].Required()) NbVerticesOnGeomVertex++;
    47964836                }
    47974837
    … …  
    48054845                nbv=0;
    48064846                for (i=0;i<Gh.nbv;i++){
    4807                         /* Add vertex only if required and not duplicate
    4808                          * (IsThe is a method of GeometricalVertex that checks
    4809                          * that we are taking into acount only one vertex is
    4810                          * 2 vertices are superimposed) */
    4811                         if (Gh[i].Required() && Gh[i].IsThe()) {//Gh  vertices Required
     4847                        /* Add vertex only if required*/
     4848                        if (Gh[i].Required()) {//Gh  vertices Required
    48124849
    48134850                                //Add the vertex (provided that nbv<maxnbv)
    … …  
    48734910                                                                else{
    48744911                                                                        e=&ei;
    4875                                                                         a=ei(0)->The();
    4876                                                                         b=ei(1)->The();
     4912                                                                        a=ei(0);
     4913                                                                        b=ei(1);
    48774914
    48784915                                                                        //check that edges has been allocated
    … …  
    49064943                                                                k=j;            // k = vertex number in edge (0 or 1)
    49074944                                                                e=&ei;          // e = reference of current edge
    4908                                                                 a=ei(k)->The(); // a = pointer toward the kth vertex of the current edge
     4945                                                                a=ei(k);        // a = pointer toward the kth vertex of the current edge
    49094946                                                                va = a->to;     // va = pointer toward vertex associated
    49104947                                                                e->SetMark();   // Mark edge
    … …  
    49144951                                                                for(;;){
    49154952                                                                        k = 1-k;
    4916                                                                         b = (*e)(k)->The();// b = pointer toward the other vertex of the current edge
     4953                                                                        b = (*e)(k);// b = pointer toward the other vertex of the current edge
    49174954                                                                        AB= b->r - a->r;   // AB = vector of the current edge
    49184955                                                                        Metric MA = background ? BTh.MetricAt(a->r) :a->m ;  //Get metric associated to A
    … …  
    54035440        }
    54045441        /*}}}1*/
    5405 /*FUNCTION Mesh::TriangleReferenceList{{{1*/
    5406 long  Mesh::TriangleReferenceList(long* reft) const {
    5407         /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ConsRefTriangle)*/
    5408 
    5409         long int verbose=0;
    5410         register Triangle *t0,*t;
    5411         register long k=0, num;   
    5412 
    5413         //initialize all triangles as -1 (outside)
    5414         for (int it=0;it<nbt;it++) reft[it]=-1;
    5415 
    5416         //loop over all subdomains
    5417         for (int i=0;i<NbSubDomains;i++){
    5418 
    5419                 //first triangle of the subdomain i
    5420                 t=t0=subdomains[i].head;
    5421 
    5422                 //check that the subdomain is not empty
    5423                 if (!t0){ ISSMERROR("At least one subdomain is empty");}
    5424 
    5425                 //loop
    5426                 do{
    5427                         k++;
    5428 
    5429                         //get current triangle number
    5430                         num = Number(t);
    5431 
    5432                         //check that num is in [0 nbt[
    5433                         if (num<0 || num>=nbt){ ISSMERROR("num<0 || num>=nbt");}
    5434 
    5435                         //reft of this triangle is the subdomain number
    5436                         reft[num]=i;
    5437 
    5438                 } while (t0 != (t=t->link));
    5439                 //stop when all triangles of subdomains have been tagged
    5440 
    5441         }
    5442         return k;   
    5443 }
    5444 /*}}}1*/
    5445 /*FUNCTION Mesh::TriangleIntNumbering{{{1*/
    5446 void Mesh::TriangleIntNumbering(long* renumbering){
    5447 
    5448         long num=0;
    5449         for (int i=0;i<nbt;i++){
    5450                 if (triangles[i].det>0) renumbering[i]=num++;
    5451                 else renumbering[i]=-1;
    5452         }
    5453         return;   
    5454 }
    5455 /*}}}1*/
    54565442
    54575443        /*Intermediary*/
    … …  
    57625748        t2->det = det2;
    57635749
    5764         t1->SetTriangleContainingTheVertex();
    5765         t2->SetTriangleContainingTheVertex();
     5750        t1->SetSingleVertexToTriangleConnectivity();
     5751        t2->SetSingleVertexToTriangleConnectivity();
    57665752} // end swap
    57675753/*}}}1*/

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

    r5120 r5124  
    107107                        void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1);
    108108                        long InsertNewPoints(long nbvold,long & NbTSwap) ;
    109                         void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
    110                         void ReNumberingVertex(long * renu);
     109                        void TrianglesRenumberBySubDomain(bool justcompress=false);
     110                        void VerticesRenumber(long * renu);
    111111                        void SmoothingVertex(int =3,double=0.3);
    112112                        Metric MetricAt (const R2 &) const;
    … …  
    120120                        long Number2(const Triangle * t) const  { return t - triangles; }
    121121                        BamgVertex* NearestVertex(Icoor1 i,Icoor1 j) ;
    122                         Triangle* FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
     122                        Triangle* TriangleFindFromCoord(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
    123123                        void ReadMesh(double* index,double* x,double* y,int nods,int nels);
    124124                        void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
    … …  
    135135
    136136                        //Inline methods
    137                         inline  void ReMakeTriangleContainingTheVertex(){
     137                        inline  void CreateSingleVertexToTriangleConnectivity(){
    138138                                for (int i=0;i<nbv;i++) vertices[i].vint=0, vertices[i].t=NULL;
    139                                 for (int i=0;i<nbt;i++) triangles[i].SetTriangleContainingTheVertex();
     139                                for (int i=0;i<nbt;i++) triangles[i].SetSingleVertexToTriangleConnectivity();
    140140                        }
    141141                        inline  void  UnMarkUnSwapTriangle(){

  • issm/trunk/src/c/objects/Bamg/Metric.h

    r3913 r5124  
    1313        class Metric;
    1414        class MatVVP2x2;
    15 
    1615
    1716        class Metric{

  • issm/trunk/src/c/objects/Bamg/R2.h

    r5120 r5124  
    77namespace bamg {
    88
    9         template <class R,class RR> class P2 {
     9        template <class R,class RR> class P2{
    1010
    1111                  public: 
    … …  
    3737          };
    3838
    39         template <class R,class RR> class P2xP2 {
     39        template <class R,class RR> class P2xP2{
    4040
    4141                  private:

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

    r5120 r5124  
    5353                        Triangle* TriangleAdj(int i) const {return TriaAdjTriangles[i&3];}
    5454                        Triangle* Quadrangle(BamgVertex * & v0,BamgVertex * & v1,BamgVertex * & v2,BamgVertex * & v3) const ;
    55                         void  ReNumbering(Triangle *tb,Triangle *te, long *renu){
     55                        void  Renumbering(Triangle *tb,Triangle *te, long *renu){
    5656                                if (link  >=tb && link  <te) link  = tb + renu[link -tb];
    5757                                if (TriaAdjTriangles[0] >=tb && TriaAdjTriangles[0] <te) TriaAdjTriangles[0] = tb + renu[TriaAdjTriangles[0]-tb];
    … …  
    5959                                if (TriaAdjTriangles[2] >=tb && TriaAdjTriangles[2] <te) TriaAdjTriangles[2] = tb + renu[TriaAdjTriangles[2]-tb];   
    6060                        }
    61                         void ReNumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
     61                        void Renumbering(BamgVertex *vb,BamgVertex *ve, long *renu){
    6262                                if (TriaVertices[0] >=vb && TriaVertices[0] <ve) TriaVertices[0] = vb + renu[TriaVertices[0]-vb];
    6363                                if (TriaVertices[1] >=vb && TriaVertices[1] <ve) TriaVertices[1] = vb + renu[TriaVertices[1]-vb];
    … …  
    8181                                }
    8282                        }
    83                         void SetTriangleContainingTheVertex() {
     83                        void SetSingleVertexToTriangleConnectivity() {
    8484                                if (TriaVertices[0]) (TriaVertices[0]->t=this,TriaVertices[0]->vint=0);
    8585                                if (TriaVertices[1]) (TriaVertices[1]->t=this,TriaVertices[1]->vint=1);

  • issm/trunk/src/c/objects/Bamg/typedefs.h

    r5120 r5124  
    1616        /*I2 and R2*/
    1717        typedef P2<Icoor1,Icoor2>  I2;
    18         typedef P2xP2<short,long>  I2xI2;
    1918        typedef P2<double,double>  R2;
    20         typedef P2<double,double>  R2xR2;
    2119}
    2220
Note: See TracChangeset for help on using the changeset viewer.