Changeset 3148

Timestamp:

03/02/10 10:14:45 (15 years ago)

Author:

Mathieu Morlighem

Message:

minor

Location:

issm/trunk/src/c/Bamgx

Files:

• Bamgx.cpp (modified) (2 diffs)
• Mesh2.h (modified) (12 diffs)
• objects/ListofIntersectionTriangles.cpp (modified) (4 diffs)
• objects/Triangle.cpp (modified) (2 diffs)
• objects/Triangles.cpp (modified) (8 diffs)

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

@@ -95 +95 @@
                 if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();
                 Th.ReNumberingTheTriangleBySubDomain();
-                //if(NbSmooth>0) Th.SmoothingVertex(NbSmooth,omega);
-                Th.SmoothingVertex(3,omega);
 
                 //Build output
@@ -159 +157 @@
                 if (verbosity>1) printf("   Generating Mesh...\n");
                 Thr=&BTh,Thb=0;
-                Triangles & Th( *(0 ?  new Triangles(*Thr,&Thr->Gh,Thb,maxnbv) :  new Triangles(maxnbv,BTh,bamgopts->KeepVertices)));
+                Triangles & Th( *(0 ?  new Triangles(*Thr,&Thr->Gh,Thb,maxnbv) :  new Triangles(maxnbv,BTh,bamgopts,bamgopts->KeepVertices)));
                 if (Thr != &BTh) delete Thr;
                 if(costheta<=1.0) Th.MakeQuadrangles(costheta);

issm/trunk/src/c/Bamgx/Mesh2.h

@@ -367 +367 @@
                 class IntersectionTriangles {
                         public: 
-                                Triangle *t;
+                                Triangle* t;
                                 Real8  bary[3];  // use if t != 0
                                 R2 x;
                                 Metric m;
-                                Real8 s;// abscisse curviline
-                                Real8 sp; // len of the previous seg in m
-                                Real8 sn;// len of the  next seg in m
+                                Real8 s; // curvilinear coordinate 
+                                Real8 sp;// length of the previous segment in m
+                                Real8 sn;// length of the next segment in m
                 };
 
@@ -380 +380 @@
                                 GeometricalEdge * e;
                                 Real8 sBegin,sEnd; // abscisse of the seg on edge parameter
-                                Real8 lBegin,lEnd; // length abscisse  set in ListofIntersectionTriangles::Length
+                                Real8 lBegin,lEnd; // length abscisse set in ListofIntersectionTriangles::Length
                                 int last;// last index  in ListofIntersectionTriangles for this Sub seg of edge
                                 R2 F(Real8 s){ 
@@ -390 +390 @@
                 };
 
-                int MaxSize; // 
-                int Size; //
-                Real8 len; //
+                int MaxSize;
+                int Size;
+                Real8 len;
                 int state;
                 IntersectionTriangles * lIntTria;
@@ -403 +403 @@
 
                 ListofIntersectionTriangles(int n=256,int m=16)
-                  :   MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
-                  NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m])  
-                {
-                 long int verbosity=0;
-                 if (verbosity>9) printf("   construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg);
+                  : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) ,
+                  NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]){
+                          long int verbosity=0;
+                          if (verbosity>9) printf("   construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg);
                 };
 
@@ -417 +416 @@
                         int NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2);
                         int NewItem(R2,const Metric & );
-                        void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1) 
-                          { 
+                        void NewSubSeg(GeometricalEdge *e,Real8 s0,Real8 s1){ 
                                 long int verbosity=0;
 
@@ -436 +434 @@
                                         lSegsI = lEn;        
                                 }
-                                if (NbSeg) 
-                                 lSegsI[NbSeg-1].last=Size;
+                                if (NbSeg) lSegsI[NbSeg-1].last=Size;
                                 lSegsI[NbSeg].e=e;
                                 lSegsI[NbSeg].sBegin=s0;
@@ -443 +440 @@
                                 NbSeg++;           
                           }
-
-                        //  void CopyMetric(int i,int j){ lIntTria[j].m=lIntTria[i].m;}
-                        //  void CopyMetric(const Metric & mm,int j){ lIntTria[j].m=mm;}
 
                         void ReShape() { 
@@ -453 +447 @@
                                         throw ErrorException(__FUNCT__,exprintf("!nw"));
                                 }
-                                for (int i=0;i<MaxSize;i++) // recopy
-                                 nw[i] = lIntTria[i];       
+                                // recopy
+                                for (int i=0;i<MaxSize;i++) nw[i] = lIntTria[i];       
                                 long int verbosity=0;
                                 if(verbosity>3) printf("   ListofIntersectionTriangles  ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg);
@@ -597 +591 @@
                 public:
 
-                        int static counter; // to kown the number of mesh in memory 
                         Geometry & Gh;   // Geometry
                         Triangles & BTh; // Background Mesh Bth==*this =>no  background 
@@ -650 +643 @@
                         Triangles(double* index,double* x,double* y,int nods,int nels);
 
-                        Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
-                                try {GeomToTriangles1(nbvx,keepBackVertices);}
+                        Triangles(Int4 nbvx,Triangles & BT,BamgOpts* bamgopts,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) {
+                                try {GeomToTriangles1(nbvx,bamgopts,keepBackVertices);}
                                 catch(...) { this->~Triangles(); throw; }
                         }
@@ -700 +693 @@
                         int SplitElement(int choice);
                         void MakeQuadTree();
-                        void NewPoints( Triangles &,int KeepVertices =1 );
+                        void NewPoints(Triangles &,BamgOpts* bamgopts,int KeepVertices=1);
                         Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; 
-                        void NewPoints(int KeepVertices=1){ NewPoints(*this,KeepVertices);}
                         void ReNumberingTheTriangleBySubDomain(bool justcompress=false);
                         void ReNumberingVertex(Int4 * renu);
@@ -744 +736 @@
                         int  CrackMesh();
                 private:
-                        void GeomToTriangles1(Int4 nbvx,int KeepVertices=1);// the real constructor mesh adaption
+                        void GeomToTriangles1(Int4 nbvx,BamgOpts* bamgopts,int KeepVertices=1);// the real constructor mesh adaption
                         void GeomToTriangles0(Int4 nbvx,BamgOpts* bamgopts);// the real constructor mesh generator
                         void PreInit(Int4);

issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp

@@ -263 +263 @@
         /*FUNCTION ListofIntersectionTriangles::Length{{{1*/
         Real8  ListofIntersectionTriangles::Length(){
+                // computation of the length
+
+                // check Size
                 if (Size<=0){
                         throw ErrorException(__FUNCT__,exprintf("Size<=0"));
                 }
-                // computation of the length      
-                R2 C;
+
                 Metric Mx,My;
                 int ii,jj;
                 R2 x,y,xy;
 
-                SegInterpolation *SegI=lSegsI;
+                SegInterpolation* SegI=lSegsI;
                 SegI=lSegsI;
-                lSegsI[NbSeg].last=Size;// improvement 
-
+                lSegsI[NbSeg].last=Size;
                 int EndSeg=Size;
 
@@ -285 +286 @@
                 for (jj=0,ii=1;ii<Size;jj=ii++) {  
                         // seg jj,ii
-                        x=y;
-                        y = lIntTria[ii].x;
+                        x  = y;
+                        y  = lIntTria[ii].x;
                         xy = y-x;
                         Mx = lIntTria[ii].m;
                         My = lIntTria[jj].m;
-                        sxy =  LengthInterpole(Mx,My,xy);
+                        sxy= LengthInterpole(Mx,My,xy);
                         s += sxy;
                         lIntTria[ii].s = s;
-                        if (ii == EndSeg) 
-                         SegI->lEnd=s,
-                                SegI++,
-                                EndSeg=SegI->last,
+                        if (ii == EndSeg){
+                                SegI->lEnd=s;
+                                SegI++;
+                                EndSeg=SegI->last;
                                 SegI->lBegin=s;
-
-                  }
+                        }
+                }
                 len = s;
                 SegI->lEnd=s;
@@ -307 +308 @@
         /*}}}1*/
         /*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/
-        Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4  nbvx){
-                long int verbosity=0;
-
-
-                const Int4 nbvold = nbv;
-                Real8 s = Length();
-                if (s <  1.5 ) return 0;
-                //////////////////////   
+        Int4 ListofIntersectionTriangles::NewPoints(Vertex* vertices,Int4 &nbv,Int4 nbvx){
+
+                //If length<1.5, do nothing
+                Real8 s=Length();
+                if (s<1.5) return 0;
+
+                const Int4 nbvold=nbv;
                 int ii = 1 ;
                 R2 y,x;
                 Metric My,Mx ;
                 Real8 sx =0,sy;
-                int nbi = Max(2,(int) (s+0.5));
-                Real8 sint = s/nbi;
-                Real8 si = sint;
+                int nbi=Max(2,(int) (s+0.5));
+                Real8 sint=s/nbi;
+                Real8 si  =sint;
 
                 int EndSeg=Size;
-                SegInterpolation *SegI=0;
-                if (NbSeg) 
-                 SegI=lSegsI,EndSeg=SegI->last;
-
-                for (int k=1;k<nbi;k++)
-                  {
-                        while ((ii < Size) && ( lIntTria[ii].s <= si )) 
-                         if (ii++ == EndSeg) 
-                          SegI++,EndSeg=SegI->last;
+                SegInterpolation* SegI=NULL;
+                if (NbSeg) SegI=lSegsI,EndSeg=SegI->last;
+
+                for (int k=1;k<nbi;k++){
+                        while ((ii < Size) && ( lIntTria[ii].s <= si )){
+                                if (ii++ == EndSeg){
+                                        SegI++;
+                                        EndSeg=SegI->last;
+                                }
+                        }
 
                         int ii1=ii-1;
@@ -347 +348 @@
                         Real8 cx = 1-cy;
                         C = SegI ? SegI->F(si): x * cx + y *cy;
+                        //C.Echo();
+                        //x.Echo();
+                        //y.Echo();
+                        //printf("cx = %g, cy=%g\n",cx,cy);
 
                         si += sint;

issm/trunk/src/c/Bamgx/objects/Triangle.cpp

@@ -98 +98 @@
         /*FUNCTION Triangle::Optim{{{1*/
         Int4  Triangle::Optim(Int2 i,int koption) {
-                // turne around in positif sens
+                // turn around (positive direction)
+                Triangle *t=this;
                 Int4 NbSwap =0;
-                Triangle  *t = this;
-                int k=0,j =OppositeEdge[i];
-                int jp = PreviousEdge[j];
-                // initialise   tp, jp the previous triangle & edge
-                Triangle *tp= TriaAdjTriangles[jp];
+                int  k = 0;
+                int  j = OppositeEdge[i];
+                int  jp= PreviousEdge[j];
+
+                // initialize tp, jp the previous triangle & edge
+                Triangle *tp=TriaAdjTriangles[jp];
                 jp = TriaAdjSharedEdge[jp]&3;
                 do {
-                        while (t->swap(j,koption))
-                          {
+                        while (t->swap(j,koption)){
+                                if (k>=20000) throw ErrorException(__FUNCT__,exprintf("k>=20000"));
                                 NbSwap++;
                                 k++;
-                                if (k>=20000){
-                                        throw ErrorException(__FUNCT__,exprintf("k>=20000"));
-                                }
                                 t=  tp->TriaAdjTriangles[jp];      // set unchange t qnd j for previous triangles
                                 j=  NextEdge[tp->TriaAdjSharedEdge[jp]&3];
-                          }
+                        }
                         // end on this  Triangle 
                         tp = t;
@@ -187 +186 @@
                                                  break;
                                          }
-                                         else 
-                                                {       
+                                         else { 
                                                  // critere de Delaunay anisotrope 
                                                  Real8 som;

issm/trunk/src/c/Bamgx/objects/Triangles.cpp

@@ -19 +19 @@
         static const  Direction NoDirOfSearch=Direction();
         long NbUnSwap =0;
-        int  Triangles::counter = 0;
 
         /*Constructors/Destructors*/
@@ -3151 +3150 @@
                 if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
                 if (verbose>3) printf("      Creating initial Constrained Delaunay Triangulation...\n");
-                if (verbose>3) printf("         Inserting points\n");
+                if (verbose>3) printf("         Inserting boundary points\n");
                 Insert();
 
@@ -3164 +3163 @@
                 // NewPointsOld(*this) ;
                 if (verbose>3) printf("      Inserting internal points\n");
-                NewPoints(*this,0) ;
+                NewPoints(*this,bamgopts,0) ;
                 if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
         }
         /*}}}1*/
         /*FUNCTION Triangles::GeomToTriangles1{{{1*/
-        void Triangles::GeomToTriangles1(Int4 inbvx,int KeepVertices) { 
+        void Triangles::GeomToTriangles1(Int4 inbvx,BamgOpts* bamgopts,int KeepVertices){ 
+
+                /*Get options*/
+                int verbosity=bamgopts->verbose;
+
                 Gh.NbRef++;// add a ref to Gh
 
@@ -3194 +3197 @@
                 }
 
-                long int verbosity=0;
                 BTh.NbRef++; // add a ref to BackGround Mesh
                 PreInit(inbvx);
@@ -3501 +3503 @@
                 FindSubDomain();
 
-                NewPoints(BTh,KeepVertices) ;
+                NewPoints(BTh,bamgopts,KeepVertices) ;
         }
         /*}}}1*/
@@ -3973 +3975 @@
 /*}}}1*/
         /*FUNCTION Triangles::NewPoints{{{1*/
-        void  Triangles::NewPoints(Triangles & Bh,int KeepVertices) {
-
-                long int verbosity=2;
-                Int4 nbtold(nbt),nbvold(nbv);
-                Int4 i,k;
-                int j;
-                Int4 *first_np_or_next_t = new Int4[nbtx];
-                Int4 NbTSwap =0;
-
-                //display info
-                if (verbosity>2)  printf("   Triangles::NewPoints\n");
-
-                /*First, insert old points*/
-
-                nbtold = nbt;
-
+        void  Triangles::NewPoints(Triangles & Bh,BamgOpts* bamgopts,int KeepVertices){
+
+                int i,j,k;
+                Int4 NbTSwap=0;
+                Int4 nbtold=nbt;
+                Int4 nbvold=nbv;
+                Int4 Headt=0;
+                Int4 next_t;
+                Int4* first_np_or_next_t=new Int4[nbtx];
+                Triangle* t=NULL;
+
+                /*Recover options*/
+                int verbosity=bamgopts->verbose;
+
+                /*First, insert old points if requested*/
                 if (KeepVertices && (&Bh != this) && (nbv+Bh.nbv< nbvx)){
-                        //   Bh.SetVertexFieldOn();
+                        if (verbosity>5) printf("         Inserting initial mesh points\n");
                         for (i=0;i<Bh.nbv;i++){ 
                                 Vertex &bv=Bh[i];
                                 if (!bv.onGeometry){
-                                        vertices[nbv].r = bv.r;
+                                        vertices[nbv].r   = bv.r;
                                         vertices[nbv++].m = bv.m;
                                 }
                         }
-                        int nbv1=nbv;
                         Bh.ReMakeTriangleContainingTheVertex();     
                         InsertNewPoints(nbvold,NbTSwap)   ;            
@@ -4004 +4004 @@
                 else Bh.ReMakeTriangleContainingTheVertex();     
 
-                Triangle *t;
                 // generation of the list of next Triangle 
-                // at 1 time we test all the triangles
-                Int4 Headt =0,next_t;
-                for(i=0;i<nbt;i++)
-                 first_np_or_next_t[i]=-(i+1);
-                // end list i >= nbt 
-                // the list of test triangle is 
-                // the next traingle on i is  -first_np_or_next_t[i]
+                for(i=0;i<nbt;i++) first_np_or_next_t[i]=-(i+1);
+                // the next traingle of i is -first_np_or_next_t[i]
+
+                // Big loop (most time consuming)
                 int iter=0;
-
-                // Big loop (most time consuming)
+                if (verbosity>5) printf("         Big loop\n");
                 do {
+                        /*Update variables*/
                         iter++;
-                        nbtold = nbt;
-                        nbvold = nbv;
-
-                        // default size of  IntersectionTriangle
+                        nbtold=nbt;
+                        nbvold=nbv;
+
+                        /*We test all triangles*/
                         i=Headt;
                         next_t=-first_np_or_next_t[i];
                         for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]){
-                                // for each triangle  t
-                                // we can change first_np_or_next_t[i]
+
+                                //check i
                                 if (i<0 || i>=nbt ){
-                                        throw ErrorException(__FUNCT__,exprintf("i<0 || i>=nbt"));
-                                }
+                                        throw ErrorException(__FUNCT__,exprintf("Index problem in NewPoints (i=%i not in [0 %i])",i,nbt-1));
+                                }
+                                //change first_np_or_next_t[i]
                                 first_np_or_next_t[i] = iter; 
+
+                                //Loop over the edges of t
                                 for(j=0;j<3;j++){
-                                        // for each edge 
                                         TriangleAdjacent tj(t,j);
-                                        Vertex & vA = * tj.EdgeVertex(0);
-                                        Vertex & vB = * tj.EdgeVertex(1);
-
-                                        if (!t->link) continue;// boundary
-                                        if (t->det <0) continue;
+                                        Vertex &vA = *tj.EdgeVertex(0);
+                                        Vertex &vB = *tj.EdgeVertex(1);
+
+                                        //if t is a boundary triangle, or tj locked, continue
+                                        if (!t->link)     continue;
+                                        if (t->det <0)    continue;
                                         if (t->Locked(j)) continue;
 
                                         TriangleAdjacent tadjj = t->Adj(j);       
-                                        Triangle * ta= tadjj;
-
-                                        if (ta->det <0) continue;         
-
-                                        R2 A = vA;
-                                        R2 B = vB;
-
+                                        Triangle* ta=tadjj;
+
+                                        //if the adjacent triangle is a boundary triangle, continur
+                                        if (ta->det<0) continue;          
+
+                                        R2 A=vA;
+                                        R2 B=vB;
                                         k=Number(ta);
 
-                                        if(first_np_or_next_t[k]==iter)  // this edge is done before 
-                                         continue; // next edge of the triangle 
-
-                                        //const Int4 NbvOld = nbv;
+                                        //if this edge has already been done, go to next edge of triangle
+                                        if(first_np_or_next_t[k]==iter) continue;
+
                                         lIntTria.SplitEdge(Bh,A,B);
                                         lIntTria.NewPoints(vertices,nbv,nbvx);
                                   } // end loop for each edge 
-
                           }// for triangle   
 
-                        if (!InsertNewPoints(nbvold,NbTSwap)) 
-                         break;
-
+                        if (!InsertNewPoints(nbvold,NbTSwap)) break;
                         for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter;
-
                         Headt = nbt; // empty list 
 
-                        // for all the triangle contening the vertex i
+                        // for all the triangle containing the vertex i
                         for (i=nbvold;i<nbv;i++){ 
                                 Vertex*          s  = vertices + i;
@@ -4085 +4079 @@
                                         ta = Next(Adj(ta));
                                 } while ( (tbegin != (Triangle*) ta)); 
-                        }   
+                        }
 
                 } while (nbv!=nbvold);