Changeset 3278


Ignore:
Timestamp:
03/15/10 14:23:55 (15 years ago)
Author:
Mathieu Morlighem
Message:

minor bug fix

Location:
issm/trunk/src/c/Bamgx
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r3267 r3278  
    5858                throw ErrorException(__FUNCT__,exprintf("maxsubdiv (%g) should be in ]1,10]",maxsubdiv));
    5959        }
     60
    6061        // no metric -> no smoothing
    6162        if (bamgopts->metric==NULL){

  • issm/trunk/src/c/Bamgx/meshtype.h

    r3277 r3278  
    118118                }
    119119        }
    120         template<class T> inline void  HeapSort(long** porder,T* c,int n){
    121                 long l,j,r,i;
    122                 T    crit;
    123                 long pos;
    124                 long* order=NULL;
    125                 order=(long*)xmalloc(n*sizeof(int));
    126                 for(i=0;i<n;i++) order[i]=i+1;
    127                 c--;                    //the array must starts at 1 and not 0
    128                 order--;
    129                 if(n<=1) return;        //return if size <=1
    130                 l=n/2+1;                //initialize l and r
    131                 r=n;
    132                 for(;;){
    133                         if(l<=1){
    134                                 crit  =c[r]; pos=order[r];
    135                                 c[r--]=c[1]; order[r+1]=order[1];
    136                                 if (r==1){
    137                                         c[1]=crit; order[1]=pos;
    138                                         order++;
    139                                         *porder=order;
    140                                         return;
    141                                 }
    142                         }
    143                         else  {crit=c[--l]; pos=order[l];}
    144                         j=l;
    145                         for(;;){
    146                                 i=j;
    147                                 j=2*j;
    148                                 if  (j>r) {c[i]=crit;order[i]=pos;break;}
    149                                 if ((j<r) && (c[j] < c[j+1]))j++;
    150                                 if (crit < c[j]) {c[i]=c[j];order[i]=order[j];}
    151                                 else{c[i]=crit;order[i]=pos;break;}
    152                         }
    153                 }
    154         }
    155120
    156121        //Inline functions

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

    r3277 r3278  
    931931        GeometricalEdge* Geometry::ProjectOnCurve(const Edge &e,double s,Vertex &V,VertexOnGeom &GV) const {
    932932                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
     933                /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/
    933934
    934935                double save_s=s;
    935936                int NbTry=0;
     937
    936938retry:
     939
    937940                s=save_s;
    938941                GeometricalEdge* on=e.onGeometry;
    939942                if (!on){
    940                         throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve error message: edge provided shonld be on geometry"));
     943                        throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve error message: edge provided should be on geometry"));
    941944                }
    942945                if (!e[0].onGeometry ||  !e[1].onGeometry){
    943946                        throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry"));
    944947                }
    945                 const Vertex &v0=e[0],&v1=e[1];
    946                 V.m = Metric(1.0-s, v0,s, v1);
    947                 const int mxe =100;
    948                 GeometricalEdge *ge[mxe+1];
    949                 int    sensge[mxe+1];
     948
     949                //Get the two vertices of the edge
     950                const Vertex &v0=e[0];
     951                const Vertex &v1=e[1];
     952
     953                //Get position of V0, V1 and vector v0->v1
     954                R2 V0=v0,V1=v1,V01=V1-V0;
     955
     956                //Get geometrical vertices corresponding to v0 and v1
     957                VertexOnGeom  vg0=*v0.onGeometry,  vg1=*v1.onGeometry;
     958
     959                //build two pointers towrad current geometrical edge
     960                GeometricalEdge *eg0=on, *eg1=on;
     961
     962                //Get edge direction and swap v0 and v1 if necessary
     963                R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
     964                int OppositeSens = (V01,AB)<0;
     965                int sens0=0,sens1=1;
     966                if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
     967
     968                //Compute metric of new vertex as a linear interpolation of the two others
     969                V.m=Metric(1.0-s,v0,s,v1);
     970
     971                const int mxe=100;
     972                GeometricalEdge* ge[mxe+1];
     973                int     sensge[mxe+1];
    950974                double  lge[mxe+1];
    951975                int bge=mxe/2,tge=bge;
    … …  
    953977                sensge[bge]=1;
    954978
    955                 R2 V0 = v0,V1=v1,V01=V1-V0;
    956                 VertexOnGeom  vg0= *v0.onGeometry,  vg1=*v1.onGeometry;
    957 
    958                 //    GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL;
    959                 GeometricalEdge * eg0=on, *eg1=on;
    960                 R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
    961                 int OppositeSens = (V01,AB) < 0;
    962                 int sens0=0,sens1=1;
    963                 if (OppositeSens)
    964                  s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
    965                 while (eg0 != (GeometricalEdge*) vg0  &&  (*eg0)(sens0) != (GeometricalVertex*) vg0){
     979                while (eg0!=(GeometricalEdge*)vg0 && (*eg0)(sens0)!=(GeometricalVertex*)vg0){
    966980                        if (bge<=0) {
    967                                 //          int kkk;
    968981                                if(NbTry) {
    969982                                        printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n");
    … …  
    975988                                  }
    976989                                NbTry++;
    977                                 goto retry;}
    978                                 GeometricalEdge* tmpge = eg0;
    979                                 ge[--bge] =eg0 = eg0->Adj[sens0];
    980                                 if (bge<0 || bge>mxe){
    981                                         throw ErrorException(__FUNCT__,exprintf("bge<0 || bge>mxe"));
    982                                 }
    983                                 sens0 = 1-( sensge[bge] = tmpge->DirAdj[sens0]);
    984                   }
     990                                goto retry;
     991                        }
     992                        GeometricalEdge* tmpge = eg0;
     993                        ge[--bge] =eg0 = eg0->Adj[sens0];
     994                        assert(bge>=0 && bge<=mxe);
     995                        sens0 = 1-( sensge[bge] = tmpge->DirAdj[sens0]);
     996                }
    985997                while (eg1 != (GeometricalEdge*) vg1  &&  (*eg1)(sens1) != (GeometricalVertex*) vg1) {
    986998                        if(tge>=mxe ) {
    … …  
    9951007                                throw ErrorException(__FUNCT__,exprintf("see above"));
    9961008                        }
    997 
    9981009                        GeometricalEdge* tmpge = eg1;
    9991010                        ge[++tge] =eg1 = eg1->Adj[sens1];
    10001011                        sensge[tge]= sens1 = 1-tmpge->DirAdj[sens1];
    1001                         if (tge<0 || tge>mxe){
    1002                                 throw ErrorException(__FUNCT__,exprintf("(tge<0 || tge>mxe)"));
    1003                         }
    1004                   }
    1005 
    1006 
    1007                 if ( (*eg0)(sens0) == (GeometricalVertex*) vg0 )
    1008                  vg0 = VertexOnGeom( *(Vertex *) vg0,*eg0,sens0);
    1009 
    1010                 if ( (*eg1)(sens1) == (GeometricalVertex*) vg1)
    1011                  vg1 = VertexOnGeom( *(Vertex *) vg1,*eg1,sens1);
     1012                        assert(tge>=0 && tge<=mxe);
     1013                }
     1014
     1015
     1016                if ((*eg0)(sens0)==(GeometricalVertex*)vg0)
     1017                 vg0=VertexOnGeom(*(Vertex*) vg0,*eg0,sens0); //vg0 = absisce
     1018
     1019                if ((*eg1)(sens1)==(GeometricalVertex*)vg1)
     1020                 vg1=VertexOnGeom(*(Vertex*) vg1,*eg1,sens1);
    10121021
    10131022                double sg;
    10141023                if (eg0 == eg1) {
    1015                         register double s0= vg0,s1=vg1;
    1016                         sg =  s0 * (1.0-s) +  s * s1;
    1017                         on=eg0;}
     1024                        register double s0=vg0,s1=vg1;
     1025                        sg =  s0*(1.0-s) +  s*s1;
     1026                        on=eg0;
     1027                }
    10181028                else {
    10191029                        R2 AA=V0,BB;
    … …  
    10221032                        double ll=0;
    10231033                        for(i=bge;i<tge;i++){
    1024                                 if ( i<0 || i>mxe){
    1025                                         throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe"));
    1026                                 }
     1034                                assert(i>=0 && i<=mxe);
    10271035                                BB =  (*ge[i])[sensge[i]];
    10281036                                lge[i]=ll += Norme2(AA-BB);
    … …  
    10431051                                                throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe"));
    10441052                                        }
    1045                                         i++,s0=1-(s1=sensge[i]),l0=l1;}
    1046                                         on=ge[i];
    1047                                         if (i==tge)
    1048                                          s1=vg1;
    1049 
    1050                                         s=(ls-l0)/(l1-l0);
    1051                                         sg =  s0 * (1.0-s) +  s * s1;   
     1053                                        i++,s0=1-(s1=sensge[i]),l0=l1;
     1054                                }
     1055                                on=ge[i];
     1056                                if (i==tge)
     1057                                 s1=vg1;
     1058
     1059                                s  =(ls-l0)/(l1-l0);
     1060                                sg =s0*(1.0-s)+s*s1;   
    10521061                }
    1053                 if (!on){
    1054                         throw ErrorException(__FUNCT__,exprintf("!on"));
    1055                 }
     1062                assert(on);
    10561063                V.r= on->F(sg);
    10571064                GV=VertexOnGeom(V,*on,sg);

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

    r3277 r3278  
    409409                                long  i1,i2;
    410410                                double s;
    411                                 i1=(long)bamgmesh->VerticesOnGeometricEdge[i*3+0]-1; //for C indexing
    412                                 i2=(long)bamgmesh->VerticesOnGeometricEdge[i*3+1]-1; //for C indexing
    413                                 s =(long)bamgmesh->VerticesOnGeometricEdge[i*3+2];
     411                                i1=(long)  bamgmesh->VerticesOnGeometricEdge[i*3+0]-1; //for C indexing
     412                                i2=(long)  bamgmesh->VerticesOnGeometricEdge[i*3+1]-1; //for C indexing
     413                                s =(double)bamgmesh->VerticesOnGeometricEdge[i*3+2];
    414414                                VerticesOnGeomEdge[i]=VertexOnGeom(vertices[i1],Gh.edges[i2],s);
    415415                        }
    … …  
    796796                                bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((Vertex*)v)+1; //back to Matlab indexing
    797797                                bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing
    798                                 bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v;
     798                                bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v; //absisce
    799799                        }
    800800                }
    … …  
    10331033                        }
    10341034                        else {
    1035                                 t->Echo();
     1035                                //t->Echo();
    10361036                                printf("\nproblem while trying to add:\n");
    10371037                                s.Echo();
    … …  
    29292929        /*}}}1*/
    29302930        /*FUNCTION Triangles::GeomToTriangles0{{{1*/
    2931         void Triangles::GeomToTriangles0(long inbvx,BamgOpts* bamgopts) {
     2931        void Triangles::GeomToTriangles0(long inbvx,BamgOpts* bamgopts){
    29322932                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/GeomToTriangles0)*/
    29332933
    … …  
    32533253
    32543254                /*Get options*/
    3255                 int verbosity=bamgopts->verbose;
     3255                int verbose=bamgopts->verbose;
    32563256
    32573257                Gh.NbRef++;// add a ref to Gh
    3258 
    32593258
    32603259                /*************************************************************************
    … …  
    327532742 internal
    32763275                 *************************************************************************/
    3277                 if (&BTh.Gh != &Gh){
    3278                         throw ErrorException(__FUNCT__,exprintf("&BTh.Gh != &Gh"));
    3279                 }
    3280 
     3276
     3277                //Check that background mesh and current mesh do have the same geometry
     3278                assert(&BTh.Gh==&Gh);
    32813279                BTh.NbRef++; // add a ref to BackGround Mesh
     3280
     3281                //Initialize new mesh
    32823282                PreInit(inbvx);
    32833283                BTh.SetVertexFieldOn();
    3284                 int * bcurve = new int[Gh.NbOfCurves]; //
     3284                int* bcurve = new int[Gh.NbOfCurves]; //
    32853285
    32863286                // we have 2 ways to make the loop
    … …  
    32953295                //  VertexOnGeom * VerticesOnGeomEdge;
    32963296
    3297 
    32983297                NbVerticesOnGeomVertex=0;
    32993298                NbVerticesOnGeomEdge=0;
    3300                 //1 copy of the  Required vertex
     3299
     3300                /*STEP 1 copy of Required vertices*/
     3301
    33013302                int i;
    3302                 for ( i=0;i<Gh.nbv;i++)
    3303                  if (Gh[i].Required()) NbVerticesOnGeomVertex++;
    3304 
    3305                 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
    3306                 VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex];
    3307                 //
    3308                 if( NbVerticesOnGeomVertex >= nbvx) {
    3309                         throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx));
    3310                 }
    3311                 if (vertices==NULL){
    3312                         throw ErrorException(__FUNCT__,exprintf("vertices==NULL"));
    3313                 }
     3303                for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
     3304                if( NbVerticesOnGeomVertex >= nbvx) { throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx));}
     3305
     3306                VerticesOnGeomVertex = new VertexOnGeom[  NbVerticesOnGeomVertex];
     3307                VertexOnBThVertex    = new VertexOnVertex[NbVerticesOnGeomVertex];
     3308
     3309                //At this point there is NO vertex but vertices should be have been allocated by PreInit
     3310                assert(vertices);
    33143311                for (i=0;i<Gh.nbv;i++){
    3315                         if (Gh[i].Required()) {//Gh  vertices Required
    3316                                 vertices[nbv] =  Gh[i];
    3317                                 vertices[nbv].i = I2(0,0);
     3312                        if (Gh[i].Required()) {//Gh vertices Required
     3313                                vertices[nbv]  =Gh[i];
     3314                                vertices[nbv].i=I2(0,0);
    33183315                                Gh[i].to = vertices + nbv;// save Geom -> Th
    33193316                                VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
    … …  
    33213318                        }
    33223319                        else Gh[i].to=0;
    3323                 }
     3320                } 
    33243321                for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){
    3325                         VertexOnGeom & vog = BTh.VerticesOnGeomVertex[i];
    3326                         if (vog.IsRequiredVertex()) {
    3327                                 GeometricalVertex * gv = vog;
     3322                        VertexOnGeom &vog=BTh.VerticesOnGeomVertex[i];
     3323                        if (vog.IsRequiredVertex()){
     3324                                GeometricalVertex* gv=vog;
    33283325                                Vertex *bv = vog;
    3329                                 if (!gv->to){// use of Geom -> Th
    3330                                         throw ErrorException(__FUNCT__,exprintf("!gv->to"));
    3331                                 }
    3332                                 VertexOnBThVertex[NbVertexOnBThVertex++] = VertexOnVertex(gv->to,bv);
     3326                                assert(gv->to); // use of Geom -> Th
     3327                                VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->to,bv);
    33333328                                gv->to->m = bv->m; // for taking the metrix of the background mesh
    3334                                 ;
    3335                         }
    3336                   }
    3337                 if (NbVertexOnBThVertex != NbVerticesOnGeomVertex){
    3338                         throw ErrorException(__FUNCT__,exprintf("NbVertexOnBThVertex != NbVerticesOnGeomVertex"));
    3339                 }
    3340                 // new stuff FH with curve
    3341                 //  find the begin of the curve in BTh
    3342                         Gh.UnMarkEdges();       
    3343                         int bfind=0;
    3344                         for (int i=0;i<Gh.NbOfCurves;i++)
    3345                           {
    3346                                 bcurve[i]=-1;
    3347                           }
    3348 
    3349                         for (int iedge=0;iedge<BTh.nbe;iedge++)
    3350                           {     
    3351                                 Edge & ei = BTh.edges[iedge];
    3352                                 for(int je=0;je<2;je++) // for the 2 extremites
    3353                                  if (!ei.onGeometry->Mark() && ei[je].onGeometry->IsRequiredVertex() )
    3354                                         {
    3355                                          // a begin of curve
    3356                                          int nc = ei.onGeometry->CurveNumber;
    3357                                          if(
    3358                                                                  ei.onGeometry==Gh.curves[nc].be    &&
    3359                                                                  (GeometricalVertex *) *ei[je].onGeometry == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] //  same extremity
    3360                                                 )     
    3361                                                 {
    3362                                                  bcurve[nc]=iedge*2+je;
    3363                                                  bfind++;       
    3364                                                 }     
     3329                        }
     3330                }
     3331                assert(NbVertexOnBThVertex==NbVerticesOnGeomVertex);
     3332
     3333                /*STEP 2: reseed bounday edges*/
     3334
     3335                //  find the begining of the curve in BTh
     3336                Gh.UnMarkEdges();       
     3337                int bfind=0;
     3338                for (int i=0;i<Gh.NbOfCurves;i++) bcurve[i]=-1;
     3339                for (int iedge=0;iedge<BTh.nbe;iedge++){     
     3340                        Edge &ei = BTh.edges[iedge];
     3341                        for(int je=0;je<2;je++){ // for the 2 extremities
     3342                                if (!ei.onGeometry->Mark() && ei[je].onGeometry->IsRequiredVertex() ){
     3343                                        // a begin of curve
     3344                                        int nc=ei.onGeometry->CurveNumber;
     3345                                        if(ei.onGeometry==Gh.curves[nc].be    &&
     3346                                                                (GeometricalVertex *) *ei[je].onGeometry == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] //  same extremity
     3347                                          ){
     3348                                                bcurve[nc]=iedge*2+je;
     3349                                                bfind++;       
    33653350                                        }
    3366                           }
    3367                         if ( bfind!=Gh.NbOfCurves){
    3368                                 throw ErrorException(__FUNCT__,exprintf("bfind!=Gh.NbOfCurves"));
    3369                         }
     3351                                }
     3352                        }
     3353                }
     3354                assert(bfind==Gh.NbOfCurves);
    33703355
    33713356                // method in 2 + 1 step
    33723357                //  0.0) compute the length and the number of vertex to do allocation
    3373                 //  1.0)  recompute the length
    3374                 //  1.1)   compute the  vertex
     3358                //  1.0) recompute the length
     3359                //  1.1) compute the  vertex
    33753360                long nbex=0,NbVerticesOnGeomEdgex=0;
    3376                 for (int step=0; step <2;step++)
    3377                   {
     3361                for (int step=0; step <2;step++){
     3362
    33783363                        long NbOfNewPoints=0;
    33793364                        long NbOfNewEdge=0;
    … …  
    33813366                        Gh.UnMarkEdges();       
    33823367                        double L=0;
    3383                         for (int icurve=0;icurve<Gh.NbOfCurves;icurve++)
    3384                           {
     3368
     3369                        for (int icurve=0;icurve<Gh.NbOfCurves;icurve++){
     3370
    33853371                                iedge=bcurve[icurve]/2;
    33863372                                int jedge=bcurve[icurve]%2;
    3387                                 if( ! Gh.curves[icurve].master) continue; // we skip all equi curve
    3388                                 Edge & ei = BTh.edges[iedge];
     3373                                if(!Gh.curves[icurve].master) continue; // we skip all equi curve
     3374                                Edge &ei = BTh.edges[iedge];
    33893375                                // warning: ei.on->Mark() can be change in
    33903376                                // loop for(jedge=0;jedge<2;jedge++)
    3391                                 // new curve 
    3392                                 // good the find a starting edge
    3393                                 double Lstep=0,Lcurve=0;// step between two points   (phase==1)
    3394                                 long NbCreatePointOnCurve=0;// Nb of new points on curve     (phase==1)
    3395 
    3396                                 for(int phase=0;phase<=step;phase++)
    3397                                   {
    3398 
    3399                                         for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next)
    3400                                           {
     3377                       
     3378                                double Lstep=0,Lcurve=0;    // step between two points   (phase==1)
     3379                                long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
     3380
     3381                                for(int phase=0;phase<=step;phase++){
     3382
     3383                                        for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next){
    34013384
    34023385                                                int icurveequi= Gh.Number(curve);
    34033386
    3404                                                 if( phase == 0 &&  icurveequi != icurve)  continue;
    3405 
    3406                                                 int k0=jedge,k1;
    3407                                                 Edge * pe=  BTh.edges+iedge;
    3408                                                 //GeometricalEdge *ong = ei.on;
    3409                                                 int iedgeequi=bcurve[icurveequi]/2;
    3410                                                 int jedgeequi=bcurve[icurveequi]%2;
     3387                                                if( phase==0 &&  icurveequi!=icurve)  continue;
     3388
     3389                                                int   k0=jedge,k1;
     3390                                                Edge* pe=  BTh.edges+iedge;
     3391                                                int   iedgeequi=bcurve[icurveequi]/2;
     3392                                                int   jedgeequi=bcurve[icurveequi]%2;
    34113393
    34123394                                                int k0equi=jedgeequi,k1equi;             
    … …  
    34143396                                                GeometricalEdge *ongequi = peequi->onGeometry;
    34153397
    3416                                                 double sNew=Lstep;// abcisse of the new points (phase==1)
     3398                                                double sNew=Lstep;// abscisse of the new points (phase==1)
    34173399                                                L=0;// length of the curve
    34183400                                                long i=0;// index of new points on the curve
    … …  
    34223404                                                Vertex *A1;
    34233405                                                VertexOnGeom *GA1;
    3424                                                 Edge * PreviousNewEdge = 0;
     3406                                                Edge* PreviousNewEdge = 0;
     3407
    34253408                                                // New Curve phase
    3426                                                 if (A0-vertices<0 || A0-vertices>=nbv){
    3427                                                         throw ErrorException(__FUNCT__,exprintf("A0-vertices<0 || A0-vertices>=nbv"));
    3428                                                 }
    3429                                                 if(ongequi->Required() )
    3430                                                   {
     3409                                                assert(A0-vertices>=0 && A0-vertices<nbv);
     3410                                                if(ongequi->Required()){
    34313411                                                        GeometricalVertex *GA1 = *(*peequi)[1-k0equi].onGeometry;
    34323412                                                        A1 = GA1->to;  //
    3433                                                   }       
    3434                                                 else
    3435                                                  for(;;){
    3436                                                          Edge &ee=*pe;
    3437                                                          Edge &eeequi=*peequi;
    3438                                                          k1 = 1-k0; // next vertex of the edge
    3439                                                          k1equi= 1 - k0equi;
    3440 
    3441                                                          if (!pe  || !ee.onGeometry){
    3442                                                                  throw ErrorException(__FUNCT__,exprintf("!pe  || !ee.on"));
    3443                                                          }
    3444                                                          ee.onGeometry->SetMark();
    3445                                                          Vertex & v0=ee[0], & v1=ee[1];
    3446                                                          R2 AB= (R2) v1 - (R2) v0;
    3447                                                          double L0=L,LAB;
    3448                                                          LAB =  LengthInterpole(v0.m,v1.m,AB);
    3449                                                          L+= LAB;   
    3450                                                          if (phase) {// computation of the new points
    3451                                                                  while ((i!=NbCreatePointOnCurve) && sNew <= L) {
    3452                                                                          if (sNew<L0){
    3453                                                                                  throw ErrorException(__FUNCT__,exprintf("sNew<L0"));
    3454                                                                          }
    3455                                                                          if (!LAB){
    3456                                                                                  throw ErrorException(__FUNCT__,exprintf("!LAB"));
    3457                                                                          }
    3458                                                                          if (!vertices || nbv>=nbvx){
    3459                                                                                  throw ErrorException(__FUNCT__,exprintf("!vertices || nbv>=nbvx"));
    3460                                                                          }
    3461                                                                          if (!edges || nbe>=nbex){
    3462                                                                                  throw ErrorException(__FUNCT__,exprintf("!edges || nbe>=nbex"));
    3463                                                                          }
    3464                                                                          if (!VerticesOnGeomEdge || NbVerticesOnGeomEdge>=NbVerticesOnGeomEdgex){
    3465                                                                                  throw ErrorException(__FUNCT__,exprintf("!VerticesOnGeomEdge || NbVerticesOnGeomEdge>=NbVerticesOnGeomEdgex"));
    3466                                                                          }
    3467                                                                          // new vertex on edge
    3468                                                                          A1=vertices+nbv++;
    3469                                                                          GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;
    3470                                                                          Edge *e = edges + nbe++;
    3471                                                                          double se= (sNew-L0)/LAB;
    3472                                                                          if (se<0 || se>=1.000000001){
    3473                                                                                  throw ErrorException(__FUNCT__,exprintf("se<0 || se>=1.000000001"));
    3474                                                                          }
    3475                                                                          se =  abscisseInterpole(v0.m,v1.m,AB,se,1);
    3476                                                                          if (se<0 || se>1){
    3477                                                                                  throw ErrorException(__FUNCT__,exprintf("se<0 || se>1"));
    3478                                                                          }
    3479                                                                          //((k1==1) != (k1==k1equi))
    3480                                                                          se = k1 ? se : 1. - se;
    3481                                                                          se = k1==k1equi ? se : 1. - se;
    3482                                                                          VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
    3483                                                                          ongequi = Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
    3484                                                                          A1->ReferenceNumber = eeequi.ref;
    3485                                                                          A1->DirOfSearch =NoDirOfSearch;
    3486                                                                          e->onGeometry = ongequi;
    3487                                                                          e->v[0]=  A0;
    3488                                                                          e->v[1]=  A1;
    3489                                                                          e->ref = eeequi.ref;
    3490                                                                          e->adj[0]=PreviousNewEdge;
    3491 
    3492                                                                          if (PreviousNewEdge)
    3493                                                                           PreviousNewEdge->adj[1] =  e;
    3494                                                                          PreviousNewEdge = e;
    3495                                                                          A0=A1;
    3496                                                                          sNew += Lstep;
    3497                                                                          if (++i== NbCreatePointOnCurve) break;
    3498                                                                  }
    3499 
    3500                                                          }               
    3501                                                          if (ee.onGeometry->CurveNumber!=ei.onGeometry->CurveNumber){
    3502                                                                  throw ErrorException(__FUNCT__,exprintf("ee.on->CurveNumber!=ei.on->CurveNumber"));
    3503                                                          }
    3504                                                          if ( ee[k1].onGeometry->IsRequiredVertex()) {
    3505                                                                  if (!eeequi[k1equi].onGeometry->IsRequiredVertex()){
    3506                                                                          throw ErrorException(__FUNCT__,exprintf("!eeequi[k1equi].on->IsRequiredVertex()"));
    3507                                                                  }
    3508                                                                  register GeometricalVertex * GA1 = *eeequi[k1equi].onGeometry;
    3509                                                                  A1=GA1->to;// the vertex in new mesh
    3510                                                                  if (A1-vertices<0 || A1-vertices>=nbv){
    3511                                                                          throw ErrorException(__FUNCT__,exprintf("A1-vertices<0 || A1-vertices>=nbv"));
    3512                                                                  }
    3513                                                                  break;}
    3514                                                                  if (!ee.adj[k1]) {
    3515                                                                          throw ErrorException(__FUNCT__,exprintf(" adj edge %i, nbe=%i, Gh.vertices=%i",BTh.Number(ee),nbe,Gh.vertices));
    3516                                                                  }
    3517                                                                  pe = ee.adj[k1]; // next edge
    3518                                                                  k0 = pe->Intersection(ee);
    3519                                                                  peequi= eeequi.adj[k1equi];  // next edge
    3520                                                                  k0equi=peequi->Intersection(eeequi);           
     3413                                                }       
     3414                                                else {
     3415                                                        for(;;){
     3416                                                                Edge &ee=*pe;
     3417                                                                Edge &eeequi=*peequi;
     3418                                                                k1 = 1-k0; // next vertex of the edge
     3419                                                                k1equi= 1 - k0equi;
     3420                                                                assert(pe && ee.onGeometry);
     3421                                                                ee.onGeometry->SetMark();
     3422                                                                Vertex & v0=ee[0], & v1=ee[1];
     3423                                                                R2 AB=(R2)v1-(R2)v0;
     3424                                                                double L0=L,LAB;
     3425                                                                LAB=LengthInterpole(v0.m,v1.m,AB);
     3426                                                                L+= LAB;
     3427
     3428                                                                if (phase){
     3429                                                                        // computation of the new points for the given curve
     3430                                                                        while ((i!=NbCreatePointOnCurve) && sNew<=L) {
     3431
     3432                                                                                //some checks
     3433                                                                                assert(sNew>=L0);
     3434                                                                                assert(LAB);
     3435                                                                                assert(vertices && nbv<nbvx);
     3436                                                                                assert(edges && nbe<nbex);
     3437                                                                                assert(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
     3438
     3439                                                                                // new vertex on edge
     3440                                                                                A1=vertices+nbv++;
     3441                                                                                GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;
     3442                                                                                Edge* e = edges + nbe++;
     3443                                                                                double se= (sNew-L0)/LAB;
     3444                                                                                if (se<0 || se>=1.000000001){
     3445                                                                                        throw ErrorException(__FUNCT__,exprintf("Problem creating point on a boundary: se=%g should be in [0 1]",se));
     3446                                                                                }
     3447                                                                                se = abscisseInterpole(v0.m,v1.m,AB,se,1);
     3448                                                                                if (se<0 || se>1){
     3449                                                                                        throw ErrorException(__FUNCT__,exprintf("Problem creating point on a boundary: se=%g should be in [0 1]",se));
     3450                                                                                }
     3451                                                                                se = k1         ? se : 1. - se;
     3452                                                                                se = k1==k1equi ? se : 1. - se;
     3453                                                                                VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
     3454                                                                                ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
     3455                                                                                A1->ReferenceNumber = eeequi.ref;
     3456                                                                                A1->DirOfSearch =NoDirOfSearch;
     3457                                                                                e->onGeometry = ongequi;
     3458                                                                                e->v[0]=A0;
     3459                                                                                e->v[1]=A1;
     3460                                                                                e->ref = eeequi.ref;
     3461                                                                                e->adj[0]=PreviousNewEdge;
     3462
     3463                                                                                if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
     3464                                                                                PreviousNewEdge=e;
     3465                                                                                A0=A1;
     3466                                                                                sNew += Lstep;
     3467                                                                                if (++i== NbCreatePointOnCurve) break;
     3468                                                                        }
     3469                                                                }
     3470
     3471                                                                //some checks
     3472                                                                assert(ee.onGeometry->CurveNumber==ei.onGeometry->CurveNumber);
     3473                                                                if (ee[k1].onGeometry->IsRequiredVertex()) {
     3474                                                                        assert(eeequi[k1equi].onGeometry->IsRequiredVertex());
     3475                                                                        register GeometricalVertex * GA1 = *eeequi[k1equi].onGeometry;
     3476                                                                        A1=GA1->to;// the vertex in new mesh
     3477                                                                        assert(A1-vertices>=0 && A1-vertices<nbv);
     3478                                                                        break;
     3479                                                                }
     3480                                                                if (!ee.adj[k1]) {
     3481                                                                        throw ErrorException(__FUNCT__,exprintf(" adj edge %i, nbe=%i, Gh.vertices=%i",BTh.Number(ee),nbe,Gh.vertices));
     3482                                                                }
     3483                                                                pe = ee.adj[k1]; // next edge
     3484                                                                k0 = pe->Intersection(ee);
     3485                                                                peequi= eeequi.adj[k1equi];  // next edge
     3486                                                                k0equi=peequi->Intersection(eeequi);           
    35213487                                                        }// for(;;) end of the curve
    3522 
    3523 
    3524                                                 if (phase) // construction of the last edge
    3525                                                   {
    3526                                                         Edge *e = edges + nbe++;
     3488                                                }
     3489
     3490
     3491                                                if (phase){ // construction of the last edge
     3492                                                        Edge* e=edges + nbe++;
    35273493                                                        e->onGeometry  = ongequi;
    3528                                                         e->v[0]=  A0;
    3529                                                         e->v[1]=        A1;
     3494                                                        e->v[0]=A0;
     3495                                                        e->v[1]=A1;
    35303496                                                        e->ref = peequi->ref;
    35313497                                                        e->adj[0]=PreviousNewEdge;
    35323498                                                        e->adj[1]=0;
    3533                                                         if (PreviousNewEdge)
    3534                                                          PreviousNewEdge->adj[1] =  e;
     3499                                                        if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
    35353500                                                        PreviousNewEdge = e;
    35363501
    3537                                                         if (i!=NbCreatePointOnCurve){
    3538                                                                 throw ErrorException(__FUNCT__,exprintf("i!=NbCreatePointOnCurve"));
    3539                                                         }
    3540 
    3541                                                   }
    3542                                           } //  end loop on equi curve
     3502                                                        assert(i==NbCreatePointOnCurve);
     3503                                                }
     3504                                        } //  end loop on equi curve
    35433505
    35443506                                        if (!phase)  { //
    … …  
    35483510                                                NbCreatePointOnCurve = NbSegOnCurve-1;
    35493511
    3550                                                 for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next)
    3551                                                   {
     3512                                                for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next){
    35523513                                                        NbOfNewEdge += NbSegOnCurve;
    35533514                                                        NbOfNewPoints += NbCreatePointOnCurve;
    3554                                                   }
    3555                                                 // do'nt
    3556                                                 //  if(NbCreatePointOnCurve<1) break;
     3515                                                }
    35573516                                        }
    3558                                   }//for(phase;;)
    3559                   } //  end of curve loop
    3560                 // end new code     
    3561                 // do the allocation
    3562                 if(step==0){
    3563                         //if(!NbOfNewPoints) break;// nothing ????? bug
    3564                         if(nbv+NbOfNewPoints > nbvx) {
    3565                                 throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,nbvx));
    3566                         }
    3567                         edges = new Edge[NbOfNewEdge];
    3568                         nbex = NbOfNewEdge;
    3569                         if(NbOfNewPoints) { //
    3570                                 VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints];
    3571                                 NbVertexOnBThEdge =NbOfNewPoints;
    3572                                 VertexOnBThEdge = new  VertexOnEdge[NbOfNewPoints];
    3573                                 NbVerticesOnGeomEdgex = NbOfNewPoints; }
     3517                                }
     3518                        }//  end of curve loop
     3519
     3520                        //Allocate memory
     3521                        if(step==0){
     3522                                if(nbv+NbOfNewPoints > nbvx) {
     3523                                        throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,nbvx));
     3524                                }
     3525                                edges = new Edge[NbOfNewEdge];
     3526                                nbex = NbOfNewEdge;
     3527                                if(NbOfNewPoints) {
     3528                                        VerticesOnGeomEdge    = new VertexOnGeom[NbOfNewPoints];
     3529                                        NbVertexOnBThEdge     = NbOfNewPoints;
     3530                                        VertexOnBThEdge       = new  VertexOnEdge[NbOfNewPoints];
     3531                                        NbVerticesOnGeomEdgex = NbOfNewPoints;
     3532                                }
    35743533                                NbOfNewPoints =0;
    35753534                                NbOfNewEdge = 0;
    3576                   }
    3577                   } // for(step;;)
    3578                 if (nbe==0){
    3579                         throw ErrorException(__FUNCT__,exprintf("nbe==0"));
    3580                 }
     3535                        }
     3536                }
     3537                assert(nbe!=0);
    35813538                delete [] bcurve;
    35823539
     3540                //Insert points inside existing triangles
     3541                if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
     3542                if (verbose>3) printf("      Creating initial Constrained Delaunay Triangulation...\n");
     3543                if (verbose>3) printf("         Inserting boundary points\n");
    35833544                Insert();
     3545
     3546                //Force the boundary
     3547                if (verbose>3) printf("         Forcing boundaries\n");
    35843548                ForceBoundary();
     3549
     3550                //Extract SubDomains
     3551                if (verbose>3) printf("         Extracting subdomains\n");
    35853552                FindSubDomain();
    35863553
     3554                if (verbose>3) printf("      Inserting internal points\n");
    35873555                NewPoints(BTh,bamgopts,KeepVertices) ;
     3556                if (verbose>4) printf("      -- current number of vertices = %i\n",nbv);
    35883557        }
    35893558        /*}}}1*/
Note: See TracChangeset for help on using the changeset viewer.