Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 3277) +++ /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 3278) @@ -58,4 +58,5 @@ throw ErrorException(__FUNCT__,exprintf("maxsubdiv (%g) should be in ]1,10]",maxsubdiv)); } + // no metric -> no smoothing if (bamgopts->metric==NULL){ Index: /issm/trunk/src/c/Bamgx/meshtype.h =================================================================== --- /issm/trunk/src/c/Bamgx/meshtype.h (revision 3277) +++ /issm/trunk/src/c/Bamgx/meshtype.h (revision 3278) @@ -118,39 +118,4 @@ } } - template inline void HeapSort(long** porder,T* c,int n){ - long l,j,r,i; - T crit; - long pos; - long* order=NULL; - order=(long*)xmalloc(n*sizeof(int)); - for(i=0;ir) {c[i]=crit;order[i]=pos;break;} - if ((j (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/ + /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/ double save_s=s; int NbTry=0; + retry: + s=save_s; GeometricalEdge* on=e.onGeometry; if (!on){ - throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve error message: edge provided shonld be on geometry")); + throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve error message: edge provided should be on geometry")); } if (!e[0].onGeometry || !e[1].onGeometry){ throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry")); } - const Vertex &v0=e[0],&v1=e[1]; - V.m = Metric(1.0-s, v0,s, v1); - const int mxe =100; - GeometricalEdge *ge[mxe+1]; - int sensge[mxe+1]; + + //Get the two vertices of the edge + const Vertex &v0=e[0]; + const Vertex &v1=e[1]; + + //Get position of V0, V1 and vector v0->v1 + R2 V0=v0,V1=v1,V01=V1-V0; + + //Get geometrical vertices corresponding to v0 and v1 + VertexOnGeom vg0=*v0.onGeometry, vg1=*v1.onGeometry; + + //build two pointers towrad current geometrical edge + GeometricalEdge *eg0=on, *eg1=on; + + //Get edge direction and swap v0 and v1 if necessary + R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; + int OppositeSens = (V01,AB)<0; + int sens0=0,sens1=1; + if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1); + + //Compute metric of new vertex as a linear interpolation of the two others + V.m=Metric(1.0-s,v0,s,v1); + + const int mxe=100; + GeometricalEdge* ge[mxe+1]; + int sensge[mxe+1]; double lge[mxe+1]; int bge=mxe/2,tge=bge; @@ -953,17 +977,6 @@ sensge[bge]=1; - R2 V0 = v0,V1=v1,V01=V1-V0; - VertexOnGeom vg0= *v0.onGeometry, vg1=*v1.onGeometry; - - // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL; - GeometricalEdge * eg0=on, *eg1=on; - R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; - int OppositeSens = (V01,AB) < 0; - int sens0=0,sens1=1; - if (OppositeSens) - s=1-s,Exchange(vg0,vg1),Exchange(V0,V1); - while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0){ + while (eg0!=(GeometricalEdge*)vg0 && (*eg0)(sens0)!=(GeometricalVertex*)vg0){ if (bge<=0) { - // int kkk; if(NbTry) { printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n"); @@ -975,12 +988,11 @@ } NbTry++; - goto retry;} - GeometricalEdge* tmpge = eg0; - ge[--bge] =eg0 = eg0->Adj[sens0]; - if (bge<0 || bge>mxe){ - throw ErrorException(__FUNCT__,exprintf("bge<0 || bge>mxe")); - } - sens0 = 1-( sensge[bge] = tmpge->DirAdj[sens0]); - } + goto retry; + } + GeometricalEdge* tmpge = eg0; + ge[--bge] =eg0 = eg0->Adj[sens0]; + assert(bge>=0 && bge<=mxe); + sens0 = 1-( sensge[bge] = tmpge->DirAdj[sens0]); + } while (eg1 != (GeometricalEdge*) vg1 && (*eg1)(sens1) != (GeometricalVertex*) vg1) { if(tge>=mxe ) { @@ -995,25 +1007,23 @@ throw ErrorException(__FUNCT__,exprintf("see above")); } - GeometricalEdge* tmpge = eg1; ge[++tge] =eg1 = eg1->Adj[sens1]; sensge[tge]= sens1 = 1-tmpge->DirAdj[sens1]; - if (tge<0 || tge>mxe){ - throw ErrorException(__FUNCT__,exprintf("(tge<0 || tge>mxe)")); - } - } - - - if ( (*eg0)(sens0) == (GeometricalVertex*) vg0 ) - vg0 = VertexOnGeom( *(Vertex *) vg0,*eg0,sens0); - - if ( (*eg1)(sens1) == (GeometricalVertex*) vg1) - vg1 = VertexOnGeom( *(Vertex *) vg1,*eg1,sens1); + assert(tge>=0 && tge<=mxe); + } + + + if ((*eg0)(sens0)==(GeometricalVertex*)vg0) + vg0=VertexOnGeom(*(Vertex*) vg0,*eg0,sens0); //vg0 = absisce + + if ((*eg1)(sens1)==(GeometricalVertex*)vg1) + vg1=VertexOnGeom(*(Vertex*) vg1,*eg1,sens1); double sg; if (eg0 == eg1) { - register double s0= vg0,s1=vg1; - sg = s0 * (1.0-s) + s * s1; - on=eg0;} + register double s0=vg0,s1=vg1; + sg = s0*(1.0-s) + s*s1; + on=eg0; + } else { R2 AA=V0,BB; @@ -1022,7 +1032,5 @@ double ll=0; for(i=bge;imxe){ - throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe")); - } + assert(i>=0 && i<=mxe); BB = (*ge[i])[sensge[i]]; lge[i]=ll += Norme2(AA-BB); @@ -1043,15 +1051,14 @@ throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe")); } - i++,s0=1-(s1=sensge[i]),l0=l1;} - on=ge[i]; - if (i==tge) - s1=vg1; - - s=(ls-l0)/(l1-l0); - sg = s0 * (1.0-s) + s * s1; + i++,s0=1-(s1=sensge[i]),l0=l1; + } + on=ge[i]; + if (i==tge) + s1=vg1; + + s =(ls-l0)/(l1-l0); + sg =s0*(1.0-s)+s*s1; } - if (!on){ - throw ErrorException(__FUNCT__,exprintf("!on")); - } + assert(on); V.r= on->F(sg); GV=VertexOnGeom(V,*on,sg); Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3277) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3278) @@ -409,7 +409,7 @@ long i1,i2; double s; - i1=(long)bamgmesh->VerticesOnGeometricEdge[i*3+0]-1; //for C indexing - i2=(long)bamgmesh->VerticesOnGeometricEdge[i*3+1]-1; //for C indexing - s =(long)bamgmesh->VerticesOnGeometricEdge[i*3+2]; + i1=(long) bamgmesh->VerticesOnGeometricEdge[i*3+0]-1; //for C indexing + i2=(long) bamgmesh->VerticesOnGeometricEdge[i*3+1]-1; //for C indexing + s =(double)bamgmesh->VerticesOnGeometricEdge[i*3+2]; VerticesOnGeomEdge[i]=VertexOnGeom(vertices[i1],Gh.edges[i2],s); } @@ -796,5 +796,5 @@ bamgmesh->VerticesOnGeometricEdge[i*3+0]=Number((Vertex*)v)+1; //back to Matlab indexing bamgmesh->VerticesOnGeometricEdge[i*3+1]=Gh.Number((const GeometricalEdge*)v)+1; //back to Matlab indexing - bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v; + bamgmesh->VerticesOnGeometricEdge[i*3+2]=(double)v; //absisce } } @@ -1033,5 +1033,5 @@ } else { - t->Echo(); + //t->Echo(); printf("\nproblem while trying to add:\n"); s.Echo(); @@ -2929,5 +2929,5 @@ /*}}}1*/ /*FUNCTION Triangles::GeomToTriangles0{{{1*/ - void Triangles::GeomToTriangles0(long inbvx,BamgOpts* bamgopts) { + void Triangles::GeomToTriangles0(long inbvx,BamgOpts* bamgopts){ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/GeomToTriangles0)*/ @@ -3253,8 +3253,7 @@ /*Get options*/ - int verbosity=bamgopts->verbose; + int verbose=bamgopts->verbose; Gh.NbRef++;// add a ref to Gh - /************************************************************************* @@ -3275,12 +3274,13 @@ 2 internal *************************************************************************/ - if (&BTh.Gh != &Gh){ - throw ErrorException(__FUNCT__,exprintf("&BTh.Gh != &Gh")); - } - + + //Check that background mesh and current mesh do have the same geometry + assert(&BTh.Gh==&Gh); BTh.NbRef++; // add a ref to BackGround Mesh + + //Initialize new mesh PreInit(inbvx); BTh.SetVertexFieldOn(); - int * bcurve = new int[Gh.NbOfCurves]; // + int* bcurve = new int[Gh.NbOfCurves]; // // we have 2 ways to make the loop @@ -3295,25 +3295,22 @@ // VertexOnGeom * VerticesOnGeomEdge; - NbVerticesOnGeomVertex=0; NbVerticesOnGeomEdge=0; - //1 copy of the Required vertex + + /*STEP 1 copy of Required vertices*/ + int i; - for ( i=0;i= nbvx) { - throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx)); - } - if (vertices==NULL){ - throw ErrorException(__FUNCT__,exprintf("vertices==NULL")); - } + for (i=0;i= nbvx) { throw ErrorException(__FUNCT__,exprintf("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 + assert(vertices); for (i=0;i Th VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]); @@ -3321,59 +3318,47 @@ } else Gh[i].to=0; - } + } for (i=0;ito){// use of Geom -> Th - throw ErrorException(__FUNCT__,exprintf("!gv->to")); - } - VertexOnBThVertex[NbVertexOnBThVertex++] = VertexOnVertex(gv->to,bv); + assert(gv->to); // use of Geom -> Th + VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->to,bv); gv->to->m = bv->m; // for taking the metrix of the background mesh - ; - } - } - if (NbVertexOnBThVertex != NbVerticesOnGeomVertex){ - throw ErrorException(__FUNCT__,exprintf("NbVertexOnBThVertex != NbVerticesOnGeomVertex")); - } - // new stuff FH with curve - // find the begin of the curve in BTh - Gh.UnMarkEdges(); - int bfind=0; - for (int i=0;iMark() && ei[je].onGeometry->IsRequiredVertex() ) - { - // a begin of curve - int nc = ei.onGeometry->CurveNumber; - if( - ei.onGeometry==Gh.curves[nc].be && - (GeometricalVertex *) *ei[je].onGeometry == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] // same extremity - ) - { - bcurve[nc]=iedge*2+je; - bfind++; - } + } + } + assert(NbVertexOnBThVertex==NbVerticesOnGeomVertex); + + /*STEP 2: reseed bounday edges*/ + + // find the begining of the curve in BTh + Gh.UnMarkEdges(); + int bfind=0; + for (int i=0;iMark() && ei[je].onGeometry->IsRequiredVertex() ){ + // a begin of curve + int nc=ei.onGeometry->CurveNumber; + if(ei.onGeometry==Gh.curves[nc].be && + (GeometricalVertex *) *ei[je].onGeometry == &(*Gh.curves[nc].be)[Gh.curves[nc].kb] // same extremity + ){ + bcurve[nc]=iedge*2+je; + bfind++; } - } - if ( bfind!=Gh.NbOfCurves){ - throw ErrorException(__FUNCT__,exprintf("bfind!=Gh.NbOfCurves")); - } + } + } + } + assert(bfind==Gh.NbOfCurves); // method in 2 + 1 step // 0.0) compute the length and the number of vertex to do allocation - // 1.0) recompute the length - // 1.1) compute the vertex + // 1.0) recompute the length + // 1.1) compute the vertex long nbex=0,NbVerticesOnGeomEdgex=0; - for (int step=0; step <2;step++) - { + for (int step=0; step <2;step++){ + long NbOfNewPoints=0; long NbOfNewEdge=0; @@ -3381,32 +3366,29 @@ Gh.UnMarkEdges(); double L=0; - for (int icurve=0;icurveMark() can be change in // loop for(jedge=0;jedge<2;jedge++) - // new curve - // good the find a starting edge - double Lstep=0,Lcurve=0;// step between two points (phase==1) - long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) - - for(int phase=0;phase<=step;phase++) - { - - for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next) - { + + double Lstep=0,Lcurve=0; // step between two points (phase==1) + long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) + + for(int phase=0;phase<=step;phase++){ + + for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next){ int icurveequi= Gh.Number(curve); - if( phase == 0 && icurveequi != icurve) continue; - - int k0=jedge,k1; - Edge * pe= BTh.edges+iedge; - //GeometricalEdge *ong = ei.on; - int iedgeequi=bcurve[icurveequi]/2; - int jedgeequi=bcurve[icurveequi]%2; + if( phase==0 && icurveequi!=icurve) continue; + + int k0=jedge,k1; + Edge* pe= BTh.edges+iedge; + int iedgeequi=bcurve[icurveequi]/2; + int jedgeequi=bcurve[icurveequi]%2; int k0equi=jedgeequi,k1equi; @@ -3414,5 +3396,5 @@ GeometricalEdge *ongequi = peequi->onGeometry; - double sNew=Lstep;// abcisse of the new points (phase==1) + double sNew=Lstep;// abscisse of the new points (phase==1) L=0;// length of the curve long i=0;// index of new points on the curve @@ -3422,123 +3404,103 @@ Vertex *A1; VertexOnGeom *GA1; - Edge * PreviousNewEdge = 0; + Edge* PreviousNewEdge = 0; + // New Curve phase - if (A0-vertices<0 || A0-vertices>=nbv){ - throw ErrorException(__FUNCT__,exprintf("A0-vertices<0 || A0-vertices>=nbv")); - } - if(ongequi->Required() ) - { + assert(A0-vertices>=0 && A0-verticesRequired()){ GeometricalVertex *GA1 = *(*peequi)[1-k0equi].onGeometry; A1 = GA1->to; // - } - else - for(;;){ - Edge &ee=*pe; - Edge &eeequi=*peequi; - k1 = 1-k0; // next vertex of the edge - k1equi= 1 - k0equi; - - if (!pe || !ee.onGeometry){ - throw ErrorException(__FUNCT__,exprintf("!pe || !ee.on")); - } - ee.onGeometry->SetMark(); - Vertex & v0=ee[0], & v1=ee[1]; - R2 AB= (R2) v1 - (R2) v0; - double L0=L,LAB; - LAB = LengthInterpole(v0.m,v1.m,AB); - L+= LAB; - if (phase) {// computation of the new points - while ((i!=NbCreatePointOnCurve) && sNew <= L) { - if (sNew=nbvx){ - throw ErrorException(__FUNCT__,exprintf("!vertices || nbv>=nbvx")); - } - if (!edges || nbe>=nbex){ - throw ErrorException(__FUNCT__,exprintf("!edges || nbe>=nbex")); - } - if (!VerticesOnGeomEdge || NbVerticesOnGeomEdge>=NbVerticesOnGeomEdgex){ - throw ErrorException(__FUNCT__,exprintf("!VerticesOnGeomEdge || NbVerticesOnGeomEdge>=NbVerticesOnGeomEdgex")); - } - // new vertex on edge - A1=vertices+nbv++; - GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge; - Edge *e = edges + nbe++; - double se= (sNew-L0)/LAB; - if (se<0 || se>=1.000000001){ - throw ErrorException(__FUNCT__,exprintf("se<0 || se>=1.000000001")); - } - se = abscisseInterpole(v0.m,v1.m,AB,se,1); - if (se<0 || se>1){ - throw ErrorException(__FUNCT__,exprintf("se<0 || se>1")); - } - //((k1==1) != (k1==k1equi)) - se = k1 ? se : 1. - se; - se = k1==k1equi ? se : 1. - se; - VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save - ongequi = Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); - A1->ReferenceNumber = eeequi.ref; - A1->DirOfSearch =NoDirOfSearch; - e->onGeometry = ongequi; - e->v[0]= A0; - e->v[1]= A1; - e->ref = eeequi.ref; - e->adj[0]=PreviousNewEdge; - - if (PreviousNewEdge) - PreviousNewEdge->adj[1] = e; - PreviousNewEdge = e; - A0=A1; - sNew += Lstep; - if (++i== NbCreatePointOnCurve) break; - } - - } - if (ee.onGeometry->CurveNumber!=ei.onGeometry->CurveNumber){ - throw ErrorException(__FUNCT__,exprintf("ee.on->CurveNumber!=ei.on->CurveNumber")); - } - if ( ee[k1].onGeometry->IsRequiredVertex()) { - if (!eeequi[k1equi].onGeometry->IsRequiredVertex()){ - throw ErrorException(__FUNCT__,exprintf("!eeequi[k1equi].on->IsRequiredVertex()")); - } - register GeometricalVertex * GA1 = *eeequi[k1equi].onGeometry; - A1=GA1->to;// the vertex in new mesh - if (A1-vertices<0 || A1-vertices>=nbv){ - throw ErrorException(__FUNCT__,exprintf("A1-vertices<0 || A1-vertices>=nbv")); - } - break;} - if (!ee.adj[k1]) { - throw ErrorException(__FUNCT__,exprintf(" adj edge %i, nbe=%i, Gh.vertices=%i",BTh.Number(ee),nbe,Gh.vertices)); - } - pe = ee.adj[k1]; // next edge - k0 = pe->Intersection(ee); - peequi= eeequi.adj[k1equi]; // next edge - k0equi=peequi->Intersection(eeequi); + } + else { + for(;;){ + Edge &ee=*pe; + Edge &eeequi=*peequi; + k1 = 1-k0; // next vertex of the edge + k1equi= 1 - k0equi; + assert(pe && ee.onGeometry); + ee.onGeometry->SetMark(); + Vertex & v0=ee[0], & v1=ee[1]; + R2 AB=(R2)v1-(R2)v0; + double L0=L,LAB; + LAB=LengthInterpole(v0.m,v1.m,AB); + L+= LAB; + + if (phase){ + // computation of the new points for the given curve + while ((i!=NbCreatePointOnCurve) && sNew<=L) { + + //some checks + assert(sNew>=L0); + assert(LAB); + assert(vertices && nbv=1.000000001){ + throw ErrorException(__FUNCT__,exprintf("Problem creating point on a boundary: se=%g should be in [0 1]",se)); + } + se = abscisseInterpole(v0.m,v1.m,AB,se,1); + if (se<0 || se>1){ + throw ErrorException(__FUNCT__,exprintf("Problem creating point on a boundary: se=%g should be in [0 1]",se)); + } + se = k1 ? se : 1. - se; + se = k1==k1equi ? se : 1. - se; + VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save + ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); + A1->ReferenceNumber = eeequi.ref; + A1->DirOfSearch =NoDirOfSearch; + e->onGeometry = ongequi; + e->v[0]=A0; + e->v[1]=A1; + e->ref = eeequi.ref; + e->adj[0]=PreviousNewEdge; + + if (PreviousNewEdge) PreviousNewEdge->adj[1]=e; + PreviousNewEdge=e; + A0=A1; + sNew += Lstep; + if (++i== NbCreatePointOnCurve) break; + } + } + + //some checks + assert(ee.onGeometry->CurveNumber==ei.onGeometry->CurveNumber); + if (ee[k1].onGeometry->IsRequiredVertex()) { + assert(eeequi[k1equi].onGeometry->IsRequiredVertex()); + register GeometricalVertex * GA1 = *eeequi[k1equi].onGeometry; + A1=GA1->to;// the vertex in new mesh + assert(A1-vertices>=0 && A1-verticesIntersection(ee); + peequi= eeequi.adj[k1equi]; // next edge + k0equi=peequi->Intersection(eeequi); }// for(;;) end of the curve - - - if (phase) // construction of the last edge - { - Edge *e = edges + nbe++; + } + + + if (phase){ // construction of the last edge + Edge* e=edges + nbe++; e->onGeometry = ongequi; - e->v[0]= A0; - e->v[1]= A1; + e->v[0]=A0; + e->v[1]=A1; e->ref = peequi->ref; e->adj[0]=PreviousNewEdge; e->adj[1]=0; - if (PreviousNewEdge) - PreviousNewEdge->adj[1] = e; + if (PreviousNewEdge) PreviousNewEdge->adj[1]=e; PreviousNewEdge = e; - if (i!=NbCreatePointOnCurve){ - throw ErrorException(__FUNCT__,exprintf("i!=NbCreatePointOnCurve")); - } - - } - } // end loop on equi curve + assert(i==NbCreatePointOnCurve); + } + } // end loop on equi curve if (!phase) { // @@ -3548,42 +3510,49 @@ NbCreatePointOnCurve = NbSegOnCurve-1; - for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next) - { + for(Curve * curve= Gh.curves+icurve;curve;curve= curve->next){ NbOfNewEdge += NbSegOnCurve; NbOfNewPoints += NbCreatePointOnCurve; - } - // do'nt - // if(NbCreatePointOnCurve<1) break; + } } - }//for(phase;;) - } // end of curve loop - // end new code - // do the allocation - if(step==0){ - //if(!NbOfNewPoints) break;// nothing ????? bug - if(nbv+NbOfNewPoints > nbvx) { - throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,nbvx)); - } - edges = new Edge[NbOfNewEdge]; - nbex = NbOfNewEdge; - if(NbOfNewPoints) { // - VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints]; - NbVertexOnBThEdge =NbOfNewPoints; - VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints]; - NbVerticesOnGeomEdgex = NbOfNewPoints; } + } + }// end of curve loop + + //Allocate memory + if(step==0){ + if(nbv+NbOfNewPoints > nbvx) { + throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,nbvx)); + } + edges = new Edge[NbOfNewEdge]; + nbex = NbOfNewEdge; + if(NbOfNewPoints) { + VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints]; + NbVertexOnBThEdge = NbOfNewPoints; + VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints]; + NbVerticesOnGeomEdgex = NbOfNewPoints; + } NbOfNewPoints =0; NbOfNewEdge = 0; - } - } // for(step;;) - if (nbe==0){ - throw ErrorException(__FUNCT__,exprintf("nbe==0")); - } + } + } + assert(nbe!=0); delete [] bcurve; + //Insert points inside existing triangles + 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 boundary points\n"); Insert(); + + //Force the boundary + if (verbose>3) printf(" Forcing boundaries\n"); ForceBoundary(); + + //Extract SubDomains + if (verbose>3) printf(" Extracting subdomains\n"); FindSubDomain(); + if (verbose>3) printf(" Inserting internal points\n"); NewPoints(BTh,bamgopts,KeepVertices) ; + if (verbose>4) printf(" -- current number of vertices = %i\n",nbv); } /*}}}1*/