Changeset 4997

Timestamp:

08/05/10 11:33:17 (15 years ago)

Author:

Mathieu Morlighem

Message:

Some improvements in bamg

Location:

issm/trunk

Files:

• src/c/objects/Bamg/Curve.h (modified) (1 diff)
• src/c/objects/Bamg/Geometry.cpp (modified) (4 diffs)
• src/c/objects/Bamg/Triangles.cpp (modified) (11 diffs)
• src/m/classes/public/bamg.m (modified) (1 diff)
• test/Validation/AnisotropicMesh/runme1_bis.m (added)
• test/Validation/AnisotropicMesh/runme2_bis.m (added)

issm/trunk/src/c/objects/Bamg/Curve.h

@@ -13 +13 @@
         class Curve {
                 public:
-                        GeometricalEdge * be,*ee; // begin et end edge
+                        GeometricalEdge *be,*ee; // begin et end edge
                         int kb,ke;  //  begin vetex and end vertex
                         Curve* next; // next curve equi to this

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

@@ -726 +726 @@
                                                                  curves[NbOfCurves].ke=k1;
                                                          }
+
                                                          GeometricalVertex *b=(*e)(k1);
                                                          if (a == b ||  b->Required() ) break;
@@ -748 +749 @@
 
                 /*clean up*/
-                delete []next_p;
-                delete []head_v;
-                delete []eangle;
+                delete [] next_p;
+                delete [] head_v;
+                delete [] eangle;
 
         }
@@ -941 +942 @@
                                 lge[tge]=ll+=Norme2(AA-V1); 
                                 // search the geometrical edge
-                                if (s>1.0){
-                                        ISSMERROR("s>1.0");
-                                }
+                                ISSMASSERT(s<=1.0);
                                 double ls= s*ll;
                                 on =0;
@@ -951 +950 @@
                                 i=bge;
                                 while (  (l1=lge[i]) < ls ) {
-                                        if (i<0 || i>mxe){
-                                                ISSMERROR("i<0 || i>mxe");
-                                        }
+                                        ISSMASSERT(i>=0 && i<=mxe);
                                         i++,s0=1-(s1=sensge[i]),l0=l1;
                                 }

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

@@ -2737 +2737 @@
                 GeometricalVertex *a,*b;
                 MeshVertex *va,*vb;
-                GeometricalEdge * e;
+                GeometricalEdge *e;
 
                 /*Get options*/
@@ -2838 +2838 @@
 
                                                                         //check that edges has been allocated
-                                                                        if (!edges){
-                                                                                ISSMERROR("edges has not been allocated...");
-                                                                        }
+                                                                        if (!edges) ISSMERROR("edges has not been allocated...");
                                                                         edges[nbe].v[0]=a->to;
                                                                         edges[nbe].v[1]=b->to;;
@@ -2854 +2852 @@
                                                 /*If Edge is not required: on a curve*/
                                                 else {
-                                                        for ( int kstep=0;kstep<=step;kstep++){
+                                                        for (int kstep=0;kstep<=step;kstep++){
                                                                 //step=0, do nothing
                                                                 //step=1, compute the length of the curve
@@ -3044 +3042 @@
 
                 /************************************************************************* 
-                // methode in 2 step
-                // 1 - compute the number of new edge to allocate
-                // 2 - construct the edge
-remark: 
-in this part we suppose to have a background mesh with the same
-geometry 
-
-To construct the discretisation of the new mesh we have to 
-rediscretize the boundary of background Mesh 
-because we have only the pointeur from the background mesh to the geometry.
-We need the abcisse of the background mesh vertices on geometry
-so a vertex is 
-0 on GeometricalVertex ;
-1 on GeometricalEdge + abcisse
-2 internal 
+                 * method in 2 steps
+                 * 1 - compute the number of new edges to allocate
+                 * 2 - construct the edges
+                 * remark:
+                 * in this part we suppose to have a background mesh with the same geometry 
+                 * 
+                 * To construct the discretization of the new mesh we have to 
+                 * rediscretize the boundary of background Mesh 
+                 * because we have only the pointeur from the background mesh to the geometry.
+                 * We need the abcisse of the background mesh vertices on geometry
+                 * so a vertex is 
+                 * 0 on GeometricalVertex ;
+                 * 1 on GeometricalEdge + abcisse
+                 * 2 internal 
                  *************************************************************************/
 
@@ -3066 +3063 @@
 
                 //Initialize new mesh
-                PreInit(inbvx);
+                this->PreInit(inbvx);
                 BTh.SetVertexFieldOn();
                 int* bcurve = new int[Gh.NbOfCurves]; // 
 
-                // we have 2 ways to make the loop 
-                // 1) on the geometry 
-                // 2) on the background mesh
-                //  if you do the loop on geometry, we don't have the pointeur on background,
-                //  and if you do the loop in background we have the pointeur on geometry
-                // so do the walk on  background
-                //  long NbVerticesOnGeomVertex;
-                //  VertexOnGeom * VerticesOnGeomVertex;  
-                //  long NbVerticesOnGeomEdge;
-                //  VertexOnGeom * VerticesOnGeomEdge;
+                /* There are 2 ways to make the loop 
+                * 1) on the geometry 
+                * 2) on the background mesh
+                *  if you do the loop on geometry, we don't have the pointeur on background,
+                *  and if you do the loop in background we have the pointeur on geometry
+                * so do the walk on  background */
 
                 NbVerticesOnGeomVertex=0;
@@ -3088 +3081 @@
                 int i; 
                 for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
-                if( NbVerticesOnGeomVertex >= nbvx) { ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx);}
+                if(NbVerticesOnGeomVertex >= nbvx) { ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx);}
 
                 VerticesOnGeomVertex = new VertexOnGeom[  NbVerticesOnGeomVertex];
                 VertexOnBThVertex    = new VertexOnVertex[NbVerticesOnGeomVertex];
 
-                //At this point there is NO vertex but vertices should be have been allocated by PreInit
+                //At this point there is NO vertex but vertices should have been allocated by PreInit
                 ISSMASSERT(vertices);
                 for (i=0;i<Gh.nbv;i++){
@@ -3117 +3110 @@
                 ISSMASSERT(NbVertexOnBThVertex==NbVerticesOnGeomVertex);
 
-                /*STEP 2: reseed bounday edges*/
+                /*STEP 2: reseed boundary edges*/
 
                 //  find the begining of the curve in BTh
@@ -3123 +3116 @@
                 int bfind=0;
                 for (int i=0;i<Gh.NbOfCurves;i++) bcurve[i]=-1; 
+
+                /*Loop over the backgrounf mesh BTh edges*/
                 for (int iedge=0;iedge<BTh.nbe;iedge++){      
                         Edge &ei = BTh.edges[iedge];
-                        for(int je=0;je<2;je++){ // for the 2 extremities
-
-                                // If one of the vertex is required we are in a new curve
+
+                        /*Loop over the 2 vertices of the current edge*/
+                        for(int je=0;je<2;je++){
+
+                                /* If one of the vertex is required we are in a new curve*/
                                 if (ei[je].onGeometry->IsRequiredVertex()){ 
 
-                                        //Get curve number
+                                        /*Get curve number*/
                                         int nc=ei.onGeometry->CurveNumber;
                                         
+                                        //printf("Dealing with curve number %i\n",nc);
+                                        //printf("edge on geometry is same as GhCurve? %s\n",(ei.onGeometry==Gh.curves[nc].be || ei.onGeometry==Gh.curves[nc].ee)?"yes":"no");
+                                        //if(ei.onGeometry==Gh.curves[nc].be || ei.onGeometry==Gh.curves[nc].ee){
+                                        //      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");
+                                        //      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");
+                                        //}
                                         //BUG FIX from original bamg
-                                        //Check that we are on the same edge and right extrimity
+                                        /*Check that we are on the same edge and right vertex (0 or 1) */
                                         if(ei.onGeometry==Gh.curves[nc].be  && (GeometricalVertex *)*ei[je].onGeometry==&(*Gh.curves[nc].be)[Gh.curves[nc].kb]){
                                                 bcurve[nc]=iedge*2+je;
@@ -3146 +3149 @@
                         }
                 } 
-                if (bfind!=Gh.NbOfCurves){
-                        ISSMERROR("problem generating number of curves (Gh.NbOfCurves=%i bfind=%i)",Gh.NbOfCurves,bfind);
-                }
+                if (bfind!=Gh.NbOfCurves) ISSMERROR("problem generating number of curves (Gh.NbOfCurves=%i bfind=%i)",Gh.NbOfCurves,bfind);
 
                 // method in 2 + 1 step 
@@ -3154 +3155 @@
                 //  1.0) recompute the length
                 //  1.1) compute the  vertex 
+
                 long nbex=0,NbVerticesOnGeomEdgex=0;
                 for (int step=0; step <2;step++){
@@ -3163 +3165 @@
                         double L=0;
 
+                        /*Go through all geometrical curve*/
                         for (int icurve=0;icurve<Gh.NbOfCurves;icurve++){ 
 
+                                /*Get edge and vertex (index) of background mesh on this curve*/
                                 iedge=bcurve[icurve]/2;
                                 int jedge=bcurve[icurve]%2;
-                                if(!Gh.curves[icurve].master) continue; // we skip all equi curve
+
+                                /*Skip if we are on a equi curve (duplicate)*/
+                                if(!Gh.curves[icurve].master) continue; 
+
+                                /*Get edge of Bth with index iedge*/
                                 Edge &ei = BTh.edges[iedge];
-                                // warning: ei.on->Mark() can be change in
-                                // loop for(jedge=0;jedge<2;jedge++) 
                         
+                                /*Initialize variables*/
                                 double Lstep=0,Lcurve=0;    // step between two points   (phase==1) 
                                 long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 
 
+                                /*Do phase 0 to step*/
                                 for(int phase=0;phase<=step;phase++){
 
-                                        for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next){
-
+                                        /*Loop over all curves from icurve till the last curve*/
+                                        for(Curve *curve= Gh.curves+icurve;curve;curve= curve->next){
+
+                                                /*Get index of current curve*/
                                                 int icurveequi= Gh.Number(curve);
 
-                                                if( phase==0 &&  icurveequi!=icurve)  continue;
+                                                /*For phase 0, check that we are at the begining of the curve only*/
+                                                if(phase==0 &&  icurveequi!=icurve)  continue;
 
                                                 int   k0=jedge,k1;

issm/trunk/src/m/classes/public/bamg.m

@@ -233 +233 @@
 
         %process geom
-        bamg_geometry=processgeometry(bamg_geometry,getfieldvalue(options,'tol',NaN),domain(1));
+        %bamg_geometry=processgeometry(bamg_geometry,getfieldvalue(options,'tol',NaN),domain(1));
 
 elseif isstruct(md.bamg),