#include #include #include #include #include "../objects.h" #include "../../include/include.h" #include "../../shared/Exceptions/exceptions.h" namespace bamg { static const Direction NoDirOfSearch=Direction(); /*Constructors/Destructors*/ /*FUNCTION Geometry::Geometry(){{{*/ Geometry::Geometry(){ Init(); } /*}}}*/ /*FUNCTION Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){{{*/ Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){ Init(); ReadGeometry(bamggeom,bamgopts); PostRead(); } /*}}}*/ /*FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{*/ Geometry::Geometry(const Geometry & Gh) { /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/Geometry)*/ long i; *this = Gh; NbRef =0; quadtree=0; vertices = nbv ? new GeomVertex[nbv] : NULL; edges = nbe ? new GeomEdge[nbe]:NULL; curves= nbcurves ? new Curve[nbcurves]:NULL; subdomains = nbsubdomains ? new GeomSubDomain[nbsubdomains]:NULL; for (i=0;i (BAMG v1.01, MeshGeom.cpp/~Geometry)*/ if(NbRef>0){ _printString_("Trying to delete geometry and NbRef>0, probably due to an error"); return;} if(vertices) delete [] vertices; vertices=0; if(edges) delete [] edges; edges=0; if(quadtree) delete quadtree; quadtree=0; if(curves) delete [] curves; curves=0;nbcurves=0; if(subdomains) delete [] subdomains;subdomains=0; Init(); } /*}}}*/ /*IO*/ /*FUNCTION Geometry::ReadGeometry{{{*/ void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){ int verbose; nbv=0; nbe=0; nbcurves=0; double Hmin = HUGE_VAL;// the infinie value int i,j,k,n,i0,i1,i2,i3; /*initialize some variables*/ verbose= bamgopts->verbose; nbv = bamggeom->VerticesSize[0]; nbe = bamggeom->EdgesSize[0]; //some checks if (bamggeom->Vertices==NULL) _error_("the domain provided does not contain any vertex"); if (bamggeom->Edges==NULL) _error_("the domain provided does not contain any edge"); //Vertices if (bamggeom->Vertices){ if(verbose>5) _printLine_(" processing Vertices"); if (bamggeom->VerticesSize[1]!=3) _error_("Vertices should have 3 columns"); vertices = new GeomVertex[nbv]; for (i=0;iVertices[i*3+0]; vertices[i].r.y=(double)bamggeom->Vertices[i*3+1]; vertices[i].ReferenceNumber=(long)bamggeom->Vertices[i*3+2]; vertices[i].DirOfSearch=NoDirOfSearch; vertices[i].color =0; vertices[i].type=0; } /*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/ pmin = vertices[0].r; pmax = vertices[0].r; for (i=0;iEdges){ R2 zerovector(0,0); double* verticeslength=NULL; if(verbose>5) _printLine_(" processing Edges"); if (bamggeom->EdgesSize[1]!=3) _error_("Edges should have 3 columns"); edges = new GeomEdge[nbe]; //initialize verticeslength (sum of the lengths of the edges holding vertex) verticeslength = new double[nbv]; for(i=0;iEdges[i*3+0]-1; //-1 for C indexing i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing edges[i].v[0]= vertices + i1; //pointer toward vertex i1 (=&vertices[i1]) edges[i].v[1]= vertices + i2; //pointer toward vertex i2 edges[i].ReferenceNumber=(long)bamggeom->Edges[i*3+2]; //get length of edge R2 x12=vertices[i2].r-vertices[i1].r; double l12=sqrt((x12,x12)); Hmin=Min(Hmin,l12); //initialize other fields edges[i].tg[0]=zerovector; edges[i].tg[1]=zerovector; edges[i].AdjVertexIndex[0] = edges[i].AdjVertexIndex[1] = -1; edges[i].Adj[0] = edges[i].Adj[1] = NULL; edges[i].type = 0; //Cracked? if (edges[i].ReferenceNumber!=1) edges[i].SetCracked(); //prepare metric vertices[i1].color++; vertices[i2].color++; verticeslength[i1] += l12; verticeslength[i2] += l12; } // definition the default of the given mesh size for (i=0;i 0) vertices[i].m=Metric(verticeslength[i] /(double) vertices[i].color); else vertices[i].m=Metric(Hmin); } delete [] verticeslength; } else{ _error_("No edges provided"); } //hVertices if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==nbv){ if(verbose>5) _printLine_(" processing hVertices"); for (i=0;i< nbv;i++){ if (!xIsNan(bamgopts->hVertices[i])){ vertices[i].m=Metric((double)bamgopts->hVertices[i]); } } } //MetricVertices if(bamgopts->metric && bamgopts->metric[0]==nbv){ if(verbose>5) _printLine_(" processing MetricVertices"); for (i=0;i< nbv;i++) { vertices[i].m = Metric((double)bamgopts->metric[i*3+0],(double)bamgopts->metric[i*3+1],(double)bamgopts->metric[i*3+2]); } } //MaxCornerAngle if (bamgopts->MaxCornerAngle){ if(verbose>5) _printLine_(" processing MaxCornerAngle"); MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180; } //TangentAtEdges if (bamggeom->TangentAtEdges){ if(verbose>5) _printString_(" processing TangentAtEdges"); if (bamggeom->TangentAtEdgesSize[1]!=4) _error_("TangentAtEdges should have 4 columns"); int n,i,j,k; R2 tg; n=bamggeom->TangentAtEdgesSize[0]; for (k=0;kTangentAtEdges[k*4+0]-1; //for C indexing j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing tg.x=bamggeom->TangentAtEdges[k*4+2]; tg.y=bamggeom->TangentAtEdges[k*4+3]; if (i<0 || i>=nbe) _error_("TangentAtEdges first index exceeds matrix dimension"); if (j!=0 && j!=1) _error_("TangentAtEdges second index should be 1 or 2 only"); edges[i].tg[j] = tg; } } //Corners if(bamggeom->Corners){ if(verbose>5) _printString_(" processing Corners"); if (bamggeom->CornersSize[1]!=1) _error_("Corners should have 1 column"); n=bamggeom->CornersSize[0]; for (i=0;iCorners[i]-1; //for C indexing if (j>nbv-1 || j<0) _error_("Bad corner definition: should in [0 " << nbv << "]"); /*Required => at the same time SetRequired and SetCorner*/ vertices[j].SetCorner(); vertices[j].SetRequired(); } } //RequiredVertices if(bamggeom->RequiredVertices){ if(verbose>5) _printLine_(" processing RequiredVertices"); if (bamggeom->RequiredVerticesSize[1]!=1) _error_("RequiredVertices should have 1 column"); n=bamggeom->RequiredVerticesSize[0]; for (i=0;iRequiredVertices[i]-1; //for C indexing if (j>nbv-1 || j<0) _error_("Bad RequiredVerticess definition: should in [0 " << nbv << "]"); vertices[j].SetRequired(); } } //RequiredEdges if(bamggeom->RequiredEdges){ if(verbose>5) _printLine_(" processing RequiredEdges"); if (bamggeom->RequiredEdgesSize[1]!=1) _error_("RequiredEdges should have 1 column"); n=bamggeom->RequiredEdgesSize[0]; for (i=0;iRequiredEdges[i]-1; //for C indexing if (j>nbe-1 || j<0) _error_("Bad RequiredEdges definition: should in [0 " << nbe << "]"); edges[j].SetRequired(); } } //SubDomain if(bamggeom->SubDomains){ if(verbose>5) _printLine_(" processing SubDomains"); if (bamggeom->SubDomainsSize[1]!=4) _error_("SubDomains should have 4 columns"); nbsubdomains=bamggeom->SubDomainsSize[0]; subdomains = new GeomSubDomain[nbsubdomains]; for (i=0;iSubDomains[i*4+0]; i1=(int)bamggeom->SubDomains[i*4+1]; i2=(int)bamggeom->SubDomains[i*4+2]; i3=(int)bamggeom->SubDomains[i*4+3]; if (i0!=2) _error_("Bad Subdomain definition: first number should be 2 (for Edges)"); if (i1>nbe || i1<=0) _error_("Bad Subdomain definition: second number should in [1 " << nbe << "] (edge number)"); subdomains[i].edge=edges + (i1-1); subdomains[i].direction = (int) i2; subdomains[i].ReferenceNumber = i3; } } } /*}}}*/ /*FUNCTION Geometry::WriteGeometry{{{*/ void Geometry::WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts){ int verbose; int nbreq=0; int nbreqv=0; int nbtan=0; int nbcracked=0; int i,count; /*Get options*/ verbose=bamgopts->verbose; /*Vertices*/ if(verbose>5) _printLine_(" writing Vertices"); bamggeom->VerticesSize[0]=nbv; bamggeom->VerticesSize[1]=3; if (nbv){ bamggeom->Vertices=xNew(3*nbv); for (i=0;iVertices[i*3+0]=vertices[i].r.x; bamggeom->Vertices[i*3+1]=vertices[i].r.y; bamggeom->Vertices[i*3+2]=vertices[i].GetReferenceNumber(); //update counters if (vertices[i].Required()) nbreqv++; } } /*Edges*/ if(verbose>5) _printLine_(" writing Edges"); bamggeom->EdgesSize[0]=nbe; bamggeom->EdgesSize[1]=3; if (nbe){ bamggeom->Edges=xNew(3*nbe); for (i=0;iEdges[i*3+0]=GetId(edges[i][0])+1; //back to Matlab indexing bamggeom->Edges[i*3+1]=GetId(edges[i][1])+1; //back to Matlab indexing bamggeom->Edges[i*3+2]=(double)edges[i].ReferenceNumber; //update counters if (edges[i].Required()) nbreq++; if (edges[i].TgA() && edges[i][0].Corner()) nbtan++; if (edges[i].TgB() && edges[i][1].Corner()) nbtan++; } } /*RequiredEdges*/ if(verbose>5) _printLine_(" writing " << nbreq << " RequiredEdges"); bamggeom->RequiredEdgesSize[0]=nbreq; bamggeom->RequiredEdgesSize[1]=1; if (nbreq){ bamggeom->RequiredEdges=xNew(1*nbreq); count=0; for (i=0;iRequiredEdges[count]=i+1; //back to Matlab indexing count=count+1; } } } //No corners /*RequiredVertices*/ if(verbose>5) _printLine_(" writing " << nbreqv << " RequiredVertices"); bamggeom->RequiredVerticesSize[0]=nbreqv; bamggeom->RequiredVerticesSize[1]=1; if (nbreqv){ bamggeom->RequiredVertices=xNew(1*nbreqv); count=0; for (i=0;iRequiredVertices[count]=i+1; //back to Matlab indexing count=count+1; } } } /*SubDomains*/ if(verbose>5) _printLine_(" writing SubDomains"); bamggeom->SubDomainsSize[0]=nbsubdomains; bamggeom->SubDomainsSize[1]=4; if (nbsubdomains){ bamggeom->SubDomains=xNew(4*nbsubdomains); for (i=0;iSubDomains[4*i+0]=2; bamggeom->SubDomains[4*i+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing bamggeom->SubDomains[4*i+2]=subdomains[i].direction; bamggeom->SubDomains[4*i+3]=subdomains[i].ReferenceNumber; } } /*TangentAtEdges*/ if(verbose>5) _printLine_(" writing TangentAtEdges"); bamggeom->TangentAtEdgesSize[0]=nbtan; bamggeom->TangentAtEdgesSize[1]=4; if (nbtan){ bamggeom->TangentAtEdges=xNew(4*nbtan); count=0; for (i=0;iTangentAtEdges[4*i+0]=i+1; //back to Matlab indexing bamggeom->TangentAtEdges[4*i+1]=1; bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[0].x; bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[0].y; } if (edges[i].TgB() && edges[i][1].Corner()){ bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing bamggeom->TangentAtEdges[4*i+1]=2; bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[1].x; bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[1].y; } count=count+1; } } } /*}}}*/ /*Methods*/ /*FUNCTION Geometry::Echo {{{*/ void Geometry::Echo(void){ _printLine_("Geometry:"); _printLine_(" nbv (number of vertices) : " << nbv); _printLine_(" nbe (number of edges) : " << nbe); _printLine_(" nbsubdomains: " << nbsubdomains); _printLine_(" nbcurves: " << nbcurves); _printLine_(" vertices: " << vertices); _printLine_(" edges: " << edges); _printLine_(" quadtree: " << quadtree); _printLine_(" subdomains: " << subdomains); _printLine_(" curves: " << curves); _printLine_(" pmin (x,y): (" << pmin.x << " " << pmin.y << ")"); _printLine_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")"); _printLine_(" coefIcoor: " << coefIcoor); _printLine_(" MaxCornerAngle: " << MaxCornerAngle); return; } /*}}}*/ /*FUNCTION Geometry::Init{{{*/ void Geometry::Init(void){ /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/ NbRef=0; nbv=0; nbe=0; quadtree=NULL; curves=NULL; edges=NULL; vertices=NULL; nbsubdomains=0; nbcurves=0; subdomains=NULL; MaxCornerAngle = 10*Pi/180; //default is 10 degres } /*}}}*/ /*FUNCTION Geometry::MinimalHmin{{{*/ double Geometry::MinimalHmin() { /* coeffIcoor = (2^30-1)/D * We cannot go beyond hmin = D/2^30 because of the quadtree * hmin is therefore approximately 2/coeffIcoor */ return 2.0/coefIcoor; }/*}}}*/ /*FUNCTION Geometry::MaximalHmax{{{*/ double Geometry::MaximalHmax() { return Max(pmax.x-pmin.x,pmax.y-pmin.y); }/*}}}*/ /*FUNCTION Geometry::GetId(const GeomVertex &t){{{*/ long Geometry::GetId(const GeomVertex & t) const { return &t - vertices; }/*}}}*/ /*FUNCTION Geometry::GetId(const GeomVertex * t){{{*/ long Geometry::GetId(const GeomVertex * t) const { return t - vertices; }/*}}}*/ /*FUNCTION Geometry::GetId(const GeomEdge & t){{{*/ long Geometry::GetId(const GeomEdge & t) const { return &t - edges; }/*}}}*/ /*FUNCTION Geometry::GetId(const GeomEdge * t){{{*/ long Geometry::GetId(const GeomEdge * t) const { return t - edges; }/*}}}*/ /*FUNCTION Geometry::GetId(const Curve * c){{{*/ long Geometry::GetId(const Curve * c) const { return c - curves; }/*}}}*/ /*FUNCTION Geometry::Containing{{{*/ GeomEdge* Geometry::Containing(const R2 P, GeomEdge * start) const { /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/Contening)*/ GeomEdge* on =start,* pon=0; // walk with the cos on geometry int counter=0; while(pon != on){ counter++; _assert_(counter<100); pon = on; R2 A= (*on)[0]; R2 B= (*on)[1]; R2 AB = B-A; R2 AP = P-A; R2 BP = P-B; if ( (AB,AP) < 0) on = on->Adj[0]; else if ( (AB,BP) > 0) on = on->Adj[1]; else return on; } return on; } /*}}}*/ /*FUNCTION Geometry::PostRead{{{*/ void Geometry::PostRead(){ /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/AfterRead)*/ long i,j,k; int jj; long *head_v = new long[nbv]; long *next_p = new long[2*nbe]; float *eangle = new float[nbe]; double eps = 1e-20; BamgQuadtree quadtree; // build quadtree to find duplicates BamgVertex *v0 = vertices; GeomVertex *v0g = (GeomVertex*) (void*)v0; k=0; //build quadtree for this geometry for (i=0;i [0 2^30-1[^2 for a quadtree of depth 30*/ vertices[i].i=R2ToI2(vertices[i].r); /*find nearest vertex already present in the quadtree (NULL if empty)*/ BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); /*if there is a vertex found that is to close to vertices[i] -> error*/ if( v && Norme1(v->r - vertices[i].r) < eps ){ _printLine_("reference numbers: " << v->ReferenceNumber << " " << vertices[i].ReferenceNumber); _printLine_("Id: " << i+1); _printLine_("Coords: ["<r.x<<" "<r.y<<"] ["<r - edges[i].v[0]->r; double lv10=Norme2(v10); //check that its length is not 0 if(lv10==0){ delete [] next_p; delete [] head_v; delete [] eangle; _error_("Length of edge " << i << " is 0"); } //compute angle in [-Pi Pi] eangle[i] = atan2(v10.y,v10.x); //build chains head_v and next_p for (j=0;j<2;j++){ long v=GetId(edges[i].v[j]); next_p[k]=head_v[v]; head_v[v]=k++; //post increment: head_v[v]=k; and then k=k+1; } } //sort head_v by order of increasing edges angle for (i=0;i=0){ ord++; long i1=n/2; // i1 = floor (n/2) -> Edge number long j1=n%2; // j1 = 1 if n is odd -> Vertex index for this edge (0 or 1) long* pn=next_p+n; //Next vertex index n=*pn; //compute angle between horizontal axis and v0->v1 float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; //exchange if the current edge angle is smaller than the previous one if (angleold > angle){ exch=1; *pn=*po; // next_p[n] = n + 1 *po=*p; // *p=n; // next_p[n+1] = n po=pn; // po now point toward pn (invert next and current) } //else, continue going to the next edge position else{ // to have : po -> p -> pn angleold=angle; // update maximum angle po=p; // po now point toward p (current position) p=pn; // p now point toward pn (next position) } } } // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i) if(ord==2) { long n1 = head_v[i]; long n2 = next_p[n1]; long i1 = n1/2, i2 = n2/2; // edge number long j1 = n1%2, j2 = n2%2; // vertex in the edge float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; float da12 = Abs(angle2-angle1); if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) { vertices[i].SetCorner() ; } // if the edge type/referencenumber a changing then is SetRequired(); if (edges[i1].type != edges[i2].type || edges[i1].Required()){ vertices[i].SetRequired(); } if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) { vertices[i].SetRequired(); } } if(ord != 2) { vertices[i].SetCorner(); } /*close the list around the vertex to have a circular loop*/ long no=-1, ne = head_v[i]; while (ne >=0) ne = next_p[no=ne]; if(no>=0) next_p[no] = head_v[i]; } /*Check that the list makes sense (we have all the time the same vertex) * and build adjacent edges*/ k =0; for (i=0;ir -edges[i].v[0]->r;// AB = vertex0 -> vertex1 double lAB=Norme2(AB); // Get length of AB double ltg2[2]={0.0}; // initialize tangent //loop over the 2 vertices of the edge for (j=0;j<2;j++) { R2 tg =edges[i].tg[j]; double ltg=Norme2(tg); //by default, tangent=[0 0] if(ltg==0){ //if the current vertex of the edge is not a corner if(!edges[i].v[j]->Corner()){ /*The tangent is set as the vector between the * previous and next vertices connected to current vertex * normed by the edge length*/ tg = edges[i].v[1-j]->r - edges[i].Adj[j]->v[1-edges[i].AdjVertexIndex[j]]->r; ltg= Norme2(tg); tg = tg *(lAB/ltg); ltg= lAB; } //else: a Corner no tangent => nothing to do } else{ //tangent has already been computed tg = tg*(lAB/ltg),ltg=lAB; } ltg2[j] = ltg; if ((tg,AB)<0) tg = -tg; edges[i].tg[j]=tg; } if (ltg2[0]!=0) edges[i].SetTgA(); if (ltg2[1]!=0) edges[i].SetTgB(); } /* generation of all curves (from corner to corner)*/ /*We proceed in 2 steps: first allocate, second build*/ for (int step=0;step<2;step++){ //unmark all edges for (i=0;iSetMark(); nb_marked_edges++; e->CurveNumber=nbcurves; GeomVertex *b=(*e)(k1); //break if we have reached the other end of the curve if (a==b || b->Required()){ if(curves){ curves[nbcurves].LastEdge=e; curves[nbcurves].LastVertexIndex=k1; } break; } //else: go to next edge (adjacent) else{ k0 = e->AdjVertexIndex[k1];// vertex in next edge e = e->Adj[k1]; // next edge } } nbcurves++; if(level) a->SetRequired(); } } } } _assert_(nb_marked_edges && nbe); //allocate if first step if(step==0) curves=new Curve[nbcurves]; } /*clean up*/ delete [] next_p; delete [] head_v; delete [] eangle; } /*}}}*/ /*FUNCTION Geometry::ProjectOnCurve {{{*/ GeomEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const { /*Original code from Frederic Hecht (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; GeomEdge* on=e.GeomEdgeHook; if (!on){ _error_("ProjectOnCurve error message: edge provided should be on geometry"); } if (!e[0].GeomEdgeHook || !e[1].GeomEdgeHook){ _error_("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry"); } //Get the two vertices of the edge const BamgVertex &v0=e[0]; const BamgVertex &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.GeomEdgeHook, vg1=*v1.GeomEdgeHook; //build two pointers towrad current geometrical edge GeomEdge *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 direction0=0,direction1=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; GeomEdge* ge[mxe+1]; int directionge[mxe+1]; double lge[mxe+1]; int bge=mxe/2,tge=bge; ge[bge] = e.GeomEdgeHook; directionge[bge]=1; while (eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){ if (bge<=0) { if(NbTry) { _printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve"); _printLine_("That bug might come from:"); _printLine_(" 1) a mesh edge containing more than " << mxe/2 << " geometrical edges"); _printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before"); _printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)"); _error_("see above"); } NbTry++; goto retry; } GeomEdge* tmpge = eg0; ge[--bge] =eg0 = eg0->Adj[direction0]; _assert_(bge>=0 && bge<=mxe); direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]); } while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) { if(tge>=mxe ) { _printLine_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)"); NbTry++; if (NbTry<2) goto retry; _printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve"); _printLine_("That bug might come from:"); _printLine_(" 1) a mesh edge contening more than " << mxe/2 << " geometrical edges"); _printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before"); _printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)"); _error_("see above"); } GeomEdge* tmpge = eg1; ge[++tge] =eg1 = eg1->Adj[direction1]; directionge[tge]= direction1 = 1-tmpge->AdjVertexIndex[direction1]; _assert_(tge>=0 && tge<=mxe); } if ((*eg0)(direction0)==(GeomVertex*)vg0) vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce if ((*eg1)(direction1)==(GeomVertex*)vg1) vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,direction1); double sg; if (eg0 == eg1) { register double s0=vg0,s1=vg1; sg = s0*(1.0-s) + s*s1; on=eg0; } else { R2 AA=V0,BB; double s0,s1; int i=bge; double ll=0; for(i=bge;i=0 && i<=mxe); BB = (*ge[i])[directionge[i]]; lge[i]=ll += Norme2(AA-BB); AA=BB ;} lge[tge]=ll+=Norme2(AA-V1); // search the geometrical edge _assert_(s<=1.0); double ls= s*ll; on =0; s0 = vg0; s1= directionge[bge]; double l0=0,l1; i=bge; while ( (l1=lge[i]) < ls ) { _assert_(i>=0 && i<=mxe); i++,s0=1-(s1=directionge[i]),l0=l1; } on=ge[i]; if (i==tge) s1=vg1; s =(ls-l0)/(l1-l0); sg =s0*(1.0-s)+s*s1; } _assert_(on); V.r= on->F(sg); GV=VertexOnGeom(V,*on,sg); return on; } /*}}}*/ /*FUNCTION Geometry::R2ToI2{{{*/ I2 Geometry::R2ToI2(const R2 & P) const { /*coefIcoor is the coefficient used for integer coordinates: * (x-pmin.x) * Icoor x = (2^30 -1) ------------ * D * where D is the longest side of the domain (direction x or y) * so that (x-pmin.x)/D is in ]0 1[ * * coefIcoor = (2^30 -1)/D */ return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) ,(Icoor1) (coefIcoor*(P.y-pmin.y)) ); }/*}}}*/ /*FUNCTION Geometry::UnMarkEdges{{{*/ void Geometry::UnMarkEdges() { for (int i=0;i