Index: /issm/trunk/src/c/objects/Bamg/Curve.h =================================================================== --- /issm/trunk/src/c/objects/Bamg/Curve.h (revision 4996) +++ /issm/trunk/src/c/objects/Bamg/Curve.h (revision 4997) @@ -13,5 +13,5 @@ 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 Index: /issm/trunk/src/c/objects/Bamg/Geometry.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/Geometry.cpp (revision 4996) +++ /issm/trunk/src/c/objects/Bamg/Geometry.cpp (revision 4997) @@ -726,4 +726,5 @@ curves[NbOfCurves].ke=k1; } + GeometricalVertex *b=(*e)(k1); if (a == b || b->Required() ) break; @@ -748,7 +749,7 @@ /*clean up*/ - delete []next_p; - delete []head_v; - delete []eangle; + delete [] next_p; + delete [] head_v; + delete [] eangle; } @@ -941,7 +942,5 @@ 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,7 +950,5 @@ 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; } Index: /issm/trunk/src/c/objects/Bamg/Triangles.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/Triangles.cpp (revision 4996) +++ /issm/trunk/src/c/objects/Bamg/Triangles.cpp (revision 4997) @@ -2737,5 +2737,5 @@ GeometricalVertex *a,*b; MeshVertex *va,*vb; - GeometricalEdge * e; + GeometricalEdge *e; /*Get options*/ @@ -2838,7 +2838,5 @@ //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,5 +2852,5 @@ /*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,19 +3042,18 @@ /************************************************************************* - // 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,18 +3063,14 @@ //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,10 +3081,10 @@ int i; for (i=0;i= 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;iIsRequiredVertex()){ - //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,7 +3149,5 @@ } } - 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,4 +3155,5 @@ // 1.0) recompute the length // 1.1) compute the vertex + long nbex=0,NbVerticesOnGeomEdgex=0; for (int step=0; step <2;step++){ @@ -3163,23 +3165,32 @@ double L=0; + /*Go through all geometrical curve*/ for (int icurve=0;icurveMark() 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; Index: /issm/trunk/src/m/classes/public/bamg.m =================================================================== --- /issm/trunk/src/m/classes/public/bamg.m (revision 4996) +++ /issm/trunk/src/m/classes/public/bamg.m (revision 4997) @@ -233,5 +233,5 @@ %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), Index: /issm/trunk/test/Validation/AnisotropicMesh/runme1_bis.m =================================================================== --- /issm/trunk/test/Validation/AnisotropicMesh/runme1_bis.m (revision 4997) +++ /issm/trunk/test/Validation/AnisotropicMesh/runme1_bis.m (revision 4997) @@ -0,0 +1,31 @@ +%test the anisotropic mesh adaptation +%function to capture = exp(-(sqrt((md.x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2; + +%create square mesh +L=1; %in m +nx=70; %numberof nodes in x direction +ny=70; +md=model; +md=squaremesh(md,L,L,nx,ny); + +%mesh adaptation loop +plotmodel(md,'data','mesh');pause(0.5); +for i=1:1 + md.vel_obs=exp(-(sqrt((md.x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2; + md.bamg=NaN; + md=bamg(md,'field',md.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',1.3,'err',10^-4); + plotmodel(md,'data','mesh');pause(0.5); +end +for i=1:1 + md.vel_obs=exp(-(sqrt((md.x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2; + md.bamg=NaN; + md=bamg(md,'field',md.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',2,'err',10^-4); + plotmodel(md,'data','mesh');pause(0.5); +end +for i=1:3, + md.vel_obs=exp(-(sqrt((md.x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2; + md.bamg=NaN; + md=bamg(md,'field',md.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',2,'err',0.005); + plotmodel(md,'data','mesh');pause(0.5); +end +md.vel_obs=exp(-(sqrt((md.x+0.1).^2+(md.y+0.1).^2)-0.75).^2*10^6)+((md.x+0.1).^2+(md.y+0.1).^2)/2; Index: /issm/trunk/test/Validation/AnisotropicMesh/runme2_bis.m =================================================================== --- /issm/trunk/test/Validation/AnisotropicMesh/runme2_bis.m (revision 4997) +++ /issm/trunk/test/Validation/AnisotropicMesh/runme2_bis.m (revision 4997) @@ -0,0 +1,43 @@ +%test the anisotropic mesh adaptation + +%create square mesh +L=1; %in m +nx=70; %numberof nodes in x direction +ny=70; +md=model; +md=squaremesh(md,L,L,nx,ny); + +%mesh adaptation loop +plotmodel(md,'data','mesh');pause(0.5); +for i=1:1 + u=4*md.x-2; v=4*md.y-2; + md.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... + +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; + md.bamg=NaN; + md=bamg(md,'field',md.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',1.3,'err',10^-4); + plotmodel(md,'data','mesh');pause(0.5); +end +for i=1:1 + u=4*md.x-2; v=4*md.y-2; + md.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... + +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; + md.bamg=NaN; + md=bamg(md,'field',md.vel_obs,'hmin',0.005,'hmax',0.3,'gradation',2,'err',10^-4); + plotmodel(md,'data','mesh');pause(0.5); +end +for i=1:3, + u=4*md.x-2; v=4*md.y-2; + md.vel_obs=tanh(30*(u.^2+v.^2-0.25)) ... + +tanh(30*((u-0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u-0.75).^2+(v+0.75).^2-0.25)) ... + +tanh(30*((u+0.75).^2+(v-0.75).^2-0.25)) +tanh(30*((u+0.75).^2+(v+0.75).^2-0.25)) ; + md.bamg=NaN; + md=bamg(md,'field',md.vel_obs,'hmin',0.001,'hmax',0.3,'gradation',2,'err',10^-3); + plotmodel(md,'data','mesh');pause(0.5); +end +md.vel_obs=tanh(30*((4*md.x-2).^2+(4*md.y-2).^2-0.25))+... + tanh(30*((4*md.x-2.75).^2+(4*md.y-1.25).^2-2.75))+... + tanh(30*((4*md.x-2.75).^2+(4*md.y-2).^2-1.25))+... + tanh(30*((4*md.x-1.25).^2+(4*md.y-2).^2-2.75))+... + tanh(30*((4*md.x-1.25).^2+(4*md.y-2).^2-1.25));