![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
source:
issm/oecreview/Archive/18296-19100/ISSM-18492-18493.diff@
19102
| Last change on this file since 19102 was 19102, checked in by , 10 years ago | |
|---|---|
| File size: 21.3 KB | |
-
../trunk-jpl/src/c/bamg/Mesh.cpp
241 241 if (Gh.NbRef>0) Gh.NbRef--; 242 242 else if (Gh.NbRef==0) delete &Gh; 243 243 } 244 if(&BTh && (&BTh != this) && false){244 if(&BTh && (&BTh != this)){ 245 245 if (BTh.NbRef>0) BTh.NbRef--; 246 246 else if (BTh.NbRef==0) delete &BTh; 247 247 } … … 4504 4504 /*Get options*/ 4505 4505 int verbose=bamgopts->verbose; 4506 4506 4507 /*Intermediaries*/ 4508 int i,k; 4509 int nbcurves = 0; 4510 int NbNewPoints,NbEdgeCurve; 4511 double lcurve,lstep,s; 4512 const int MaxSubEdge = 10; 4507 Gh.NbRef++;// add a ref to Gh 4513 4508 4514 R2 AB; 4515 GeomVertex *a, *b; 4516 BamgVertex *va,*vb; 4517 GeomEdge *e; 4509 /************************************************************************* 4510 * method in 2 steps 4511 * 1 - compute the number of new edges to allocate 4512 * 2 - construct the edges 4513 * remark: 4514 * in this part we suppose to have a background mesh with the same geometry 4515 * 4516 * To construct the discretization of the new mesh we have to 4517 * rediscretize the boundary of background Mesh 4518 * because we have only the pointeur from the background mesh to the geometry. 4519 * We need the abcisse of the background mesh vertices on geometry 4520 * so a vertex is 4521 * 0 on GeomVertex ; 4522 * 1 on GeomEdge + abcisse 4523 * 2 internal 4524 *************************************************************************/ 4518 4525 4519 // add a ref to GH to make sure that it is not destroyed by mistake 4520 Gh.NbRef++; 4526 //Check that background mesh and current mesh do have the same geometry 4527 _assert_(&BTh.Gh==&Gh); 4528 BTh.NbRef++; // add a ref to BackGround Mesh 4521 4529 4522 //build background mesh flag (1 if background, else 0) 4523 bool background=(&BTh != this); 4530 //Initialize new mesh 4531 BTh.SetVertexFieldOn(); 4532 int* bcurve = new int[Gh.nbcurves]; // 4524 4533 4525 /*Build VerticesOnGeomVertex*/ 4534 /* There are 2 ways to make the loop 4535 * 1) on the geometry 4536 * 2) on the background mesh 4537 * if you do the loop on geometry, we don't have the pointeur on background, 4538 * and if you do the loop in background we have the pointeur on geometry 4539 * so do the walk on background */ 4526 4540 4527 //Compute the number of geometrical vertices that we are going to use to mesh 4528 for (i=0;i<Gh.nbv;i++){ 4529 if (Gh[i].Required()) NbVerticesOnGeomVertex++; 4530 } 4531 //allocate 4532 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 4533 if(NbVerticesOnGeomVertex >= maxnbv) _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv); 4534 _assert_(nbv==0); 4535 //Build VerticesOnGeomVertex 4536 for (i=0;i<Gh.nbv;i++){ 4537 /* Add vertex only if required*/ 4538 if (Gh[i].Required()) {//Gh vertices Required 4541 NbVerticesOnGeomVertex=0; 4542 NbVerticesOnGeomEdge=0; 4539 4543 4540 //Add the vertex 4541 _assert_(nbv<maxnbv); 4542 vertices[nbv]=Gh[i]; 4544 /*STEP 1 copy of Required vertices*/ 4543 4545 4544 //Add pointer from geometry (Gh) to vertex from mesh (Th) 4545 Gh[i].MeshVertexHook=vertices+nbv; 4546 int i; 4547 for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++; 4548 printf("\n"); 4549 if(NbVerticesOnGeomVertex >= maxnbv){ 4550 _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv); 4551 } 4546 4552 4547 //Build VerticesOnGeomVertex for current point4548 VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].MeshVertexHook,Gh[i]);4553 VerticesOnGeomVertex = new VertexOnGeom[ NbVerticesOnGeomVertex]; 4554 VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex]; 4549 4555 4550 //nbv increment 4556 //At this point there is NO vertex but vertices should have been allocated by Init 4557 _assert_(vertices); 4558 for (i=0;i<Gh.nbv;i++){ 4559 if (Gh[i].Required()) {//Gh vertices Required 4560 vertices[nbv] =Gh[i]; 4561 vertices[nbv].i=I2(0,0); 4562 Gh[i].MeshVertexHook = vertices + nbv;// save Geom -> Th 4563 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]); 4551 4564 nbv++; 4552 4565 } 4566 else Gh[i].MeshVertexHook=0; 4567 } 4568 for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){ 4569 VertexOnGeom &vog=BTh.VerticesOnGeomVertex[i]; 4570 if (vog.IsRequiredVertex()){ 4571 GeomVertex* gv=vog; 4572 BamgVertex *bv = vog; 4573 _assert_(gv->MeshVertexHook); // use of Geom -> Th 4574 VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->MeshVertexHook,bv); 4575 gv->MeshVertexHook->m = bv->m; // for taking the metric of the background mesh 4576 } 4553 4577 } 4578 _assert_(NbVertexOnBThVertex==NbVerticesOnGeomVertex); /*This might be due to MaxCornerAngle too small*/ 4554 4579 4555 /* Build VerticesOnGeomEdge*/4580 /*STEP 2: reseed boundary edges*/ 4556 4581 4557 //check that edges is still empty (Init) 4558 _assert_(!edges); 4582 // find the begining of the curve in BTh 4583 Gh.UnMarkEdges(); 4584 int bfind=0; 4585 for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1; 4559 4586 4560 /* Now we are going to create the first edges corresponding 4561 * to the one present in the geometry provided. 4562 * We proceed in 2 steps 4563 * -step 0: we count all the edges 4564 * we allocate the number of edges at the end of step 0 4565 * -step 1: the edges are created */ 4566 for (int step=0;step<2;step++){ 4587 /*Loop over the backgrounf mesh BTh edges*/ 4588 for (int iedge=0;iedge<BTh.nbe;iedge++){ 4589 Edge &ei = BTh.edges[iedge]; 4567 4590 4568 //initialize number of edges and number of edges max 4569 long nbex=0; 4570 nbe=0; 4571 long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge; 4572 Gh.UnMarkEdges(); 4573 nbcurves=0; 4591 /*Loop over the 2 vertices of the current edge*/ 4592 for(int je=0;je<2;je++){ 4574 4593 4575 //go through the edges of the geometry4576 for (i=0;i<Gh.nbe;i++){4594 /* If one of the vertex is required we are in a new curve*/ 4595 if (ei[je].GeomEdgeHook->IsRequiredVertex()){ 4577 4596 4578 //ei = current Geometrical edge4579 GeomEdge &ei=Gh.edges[i];4597 /*Get curve number*/ 4598 int nc=ei.GeomEdgeHook->CurveNumber; 4580 4599 4581 //loop over the two vertices of the edge ei 4582 for(int j=0;j<2;j++) { 4600 //_printf_("Dealing with curve number " << nc << "\n"); 4601 //_printf_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no\n"); 4602 //if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge){ 4603 // _printf_("Do we have the right extremity? curve first vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex])?"yes":"no\n"); 4604 // _printf_("Do we have the right extremity? curve last vertex -> " << ((GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex])?"yes":"no\n"); 4605 //} 4606 //BUG FIX from original bamg 4607 /*Check that we are on the same edge and right vertex (0 or 1) */ 4608 if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex]){ 4609 bcurve[nc]=iedge*2+je; 4610 bfind++; 4611 } 4612 else if ((ei.GeomEdgeHook==Gh.curves[nc].LastEdge && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex]) && bcurve[nc]==-1){ 4613 bcurve[nc]=iedge*2+je; 4614 bfind++; 4615 } 4616 } 4617 } 4618 } 4619 if (bfind!=Gh.nbcurves){ 4620 delete [] bcurve; 4621 _error_("problem generating number of curves (" << Gh.nbcurves << " found in the geometry but " << bfind << " curve found in the mesh)"); 4622 } 4583 4623 4584 /*Take only required vertices (corner->beginning of a new curve)*/ 4585 if (!ei.Mark() && ei[j].Required()){ 4624 // method in 2 + 1 step 4625 // 0.0) compute the length and the number of vertex to do allocation 4626 // 1.0) recompute the length 4627 // 1.1) compute the vertex 4586 4628 4587 long nbvend=0; 4588 Edge* PreviousNewEdge=NULL; 4589 lstep = -1; 4629 long nbex=0,NbVerticesOnGeomEdgex=0; 4630 for (int step=0; step <2;step++){ 4590 4631 4591 /*If Edge is required (do that only once for the 2 vertices)*/ 4592 if(ei.Required()){ 4593 if (j==0){ 4594 //do not create internal points if required (take it as is) 4595 if(step==0) nbe++; 4596 else{ 4597 e=&ei; 4598 a=ei(0); 4599 b=ei(1); 4632 long NbOfNewPoints=0; 4633 long NbOfNewEdge=0; 4634 long iedge; 4635 Gh.UnMarkEdges(); 4636 double L=0; 4600 4637 4601 //check that edges has been allocated 4602 _assert_(edges); 4603 edges[nbe].v[0]=a->MeshVertexHook; 4604 edges[nbe].v[1]=b->MeshVertexHook;; 4605 edges[nbe].ReferenceNumber = e->ReferenceNumber; 4606 edges[nbe].GeomEdgeHook = e; 4607 edges[nbe].adj[0] = 0; 4608 edges[nbe].adj[1] = 0; 4609 nbe++; 4610 } 4611 } 4612 } 4638 /*Go through all geometrical curve*/ 4639 for (int icurve=0;icurve<Gh.nbcurves;icurve++){ 4613 4640 4614 /*If Edge is not required: we are on a curve*/ 4615 else { 4616 for (int kstep=0;kstep<=step;kstep++){ 4617 //kstep=0, compute number of edges (discretize curve) 4618 //kstep=1 create the points and edge 4619 PreviousNewEdge=0; 4620 NbNewPoints=0; 4621 NbEdgeCurve=0; 4622 if (nbvend>=maxnbv) _error_("maximum number of vertices too low! Check the domain outline or increase maxnbv"); 4623 lcurve =0; 4624 s = lstep; //-1 initially, then length of each sub edge 4641 /*Get edge and vertex (index) of background mesh on this curve*/ 4642 iedge=bcurve[icurve]/2; 4643 int jedge=bcurve[icurve]%2; 4625 4644 4626 /*reminder: i = edge number, j=[0;1] vertex index in edge*/ 4627 k=j; // k = vertex index in edge (0 or 1) 4628 e=&ei; // e = reference of current edge 4629 a=ei(k); // a = pointer toward the kth vertex of the current edge 4630 va = a->MeshVertexHook; // va = pointer toward mesh vertex associated 4631 e->SetMark(); // Mark edge 4645 /*Get edge of Bth with index iedge*/ 4646 Edge &ei = BTh.edges[iedge]; 4632 4647 4633 /*Loop until we reach the end of the curve*/ 4634 for(;;){ 4635 k = 1-k; // other vertx index of the curve 4636 b = (*e)(k); // b = pointer toward the other vertex of the current edge 4637 AB= b->r - a->r; // AB = vector of the current edge 4638 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A 4639 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B 4640 double ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric 4648 /*Initialize variables*/ 4649 double Lstep=0; // step between two points (phase==1) 4650 long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 4641 4651 4642 /* We are now creating the mesh edges from the geometrical edge selected above. 4643 * The edge will be divided according to the metric previously computed and cannot 4644 * be divided more than 10 times (MaxSubEdge). */ 4652 /*Do phase 0 to step*/ 4653 for(int phase=0;phase<=step;phase++){ 4645 4654 4646 //By default, there is only one subedge that is the geometrical edge itself4647 int NbSubEdge = 1;4655 /*Current curve pointer*/ 4656 Curve *curve= Gh.curves+icurve; 4648 4657 4649 //initialize lSubEdge, holding the length of each subedge (cannot be higher than 10)4650 double lSubEdge[MaxSubEdge];4658 /*Get index of current curve*/ 4659 int icurveequi= Gh.GetId(curve); 4651 4660 4652 //Build Subedges according to the edge length 4653 if (ledge < 1.5){ 4654 //if ledge < 1.5 (between one and 2), take the edge as is 4655 lSubEdge[0] = ledge; 4656 } 4657 //else, divide the edge 4658 else { 4659 //compute number of subedges (division of the edge), Maximum is 10 4660 NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5)); 4661 /*Now, we are going to divide the edge according to the metric. 4662 * Get segment by sement along the edge. 4663 * Build lSubEdge, which holds the distance between the first vertex 4664 * of the edge and the next point on the edge according to the 4665 * discretization (each SubEdge is AB)*/ 4666 R2 A,B; 4667 A=a->r; 4668 Metric MAs=MA,MBs; 4669 ledge=0; 4670 double x =0, xstep= 1./NbSubEdge; 4671 for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { 4672 x += xstep; 4673 B = e->F(k ? x : 1-x); 4674 MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB); 4675 AB = A-B; 4676 lSubEdge[kk]=(ledge+=(MAs(AB)+MBs(AB))/2); 4677 } 4678 } 4661 /*For phase 0, check that we are at the begining of the curve only*/ 4662 if(phase==0 && icurveequi!=icurve) continue; 4679 4663 4680 double lcurveb = lcurve+ledge; 4664 int k0=jedge,k1; 4665 Edge* pe= BTh.edges+iedge; 4666 int iedgeequi=bcurve[icurveequi]/2; 4667 int jedgeequi=bcurve[icurveequi]%2; 4681 4668 4682 /*Now, create corresponding points*/ 4683 while(s>=lcurve && s<=lcurveb && nbv<nbvend){ 4669 int k0equi=jedgeequi,k1equi; 4670 Edge * peequi= BTh.edges+iedgeequi; 4671 GeomEdge *ongequi = peequi->GeomEdgeHook; 4684 4672 4685 /*Schematic of current curve 4686 * 4687 * a vb b // vertex 4688 * 0 ll0 ll1 ledge // length from a 4689 * + --- + - ... - + --S-- + --- + - ... - + // where is S 4690 * 0 kk0 kk1 NbSubEdge // Sub edge index 4691 * 4692 */ 4673 double sNew=Lstep;// abscisse of the new points (phase==1) 4674 L=0;// length of the curve 4675 long i=0;// index of new points on the curve 4676 GeomVertex * GA0 = *(*peequi)[k0equi].GeomEdgeHook; 4677 BamgVertex *A0; 4678 A0 = GA0->MeshVertexHook; // the vertex in new mesh 4679 BamgVertex *A1; 4680 VertexOnGeom *GA1; 4681 Edge* PreviousNewEdge = 0; 4693 4682 4694 double ss = s-lcurve; 4683 // New Curve phase 4684 _assert_(A0-vertices>=0 && A0-vertices<nbv); 4685 if(ongequi->Required()){ 4686 GeomVertex *GA1 = *(*peequi)[1-k0equi].GeomEdgeHook; 4687 A1 = GA1->MeshVertexHook; // 4688 } 4689 else { 4690 for(;;){ 4691 Edge &ee=*pe; 4692 Edge &eeequi=*peequi; 4693 k1 = 1-k0; // next vertex of the edge 4694 k1equi= 1 - k0equi; 4695 _assert_(pe && ee.GeomEdgeHook); 4696 ee.GeomEdgeHook->SetMark(); 4697 BamgVertex & v0=ee[0], & v1=ee[1]; 4698 R2 AB=(R2)v1-(R2)v0; 4699 double L0=L,LAB; 4700 LAB=LengthInterpole(v0.m,v1.m,AB); 4701 L+= LAB; 4695 4702 4696 /*Find the SubEdge containing ss using Dichotomy*/ 4697 int kk0=-1,kk1=NbSubEdge-1,kkk; 4698 double ll0=0,ll1=ledge,llk; 4699 while (kk1-kk0>1){ 4700 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) 4701 kk1=kkk,ll1=llk; 4702 else 4703 kk0=kkk,ll0=llk; 4704 } 4705 _assert_(kk1!=kk0); 4703 if (phase){ 4704 // computation of the new points for the given curve 4705 while ((i!=NbCreatePointOnCurve) && sNew<=L) { 4706 4706 4707 /*Curvilinear coordinate in [0 1] of ss in current edge*/ 4708 // WARNING: This is what we would do 4709 // ssa = (ss-ll0)/(ll1-ll0); 4710 // aa = (kk0+ssa)/NbSubEdge 4711 // This is what Bamg does: 4712 double sbb = (ss-ll0)/(ll1-ll0); 4713 /*Curvilinear coordinate in [0 1] of ss in current curve*/ 4714 double bb = (kk1+sbb)/NbSubEdge; 4715 double aa = 1-bb; 4707 //some checks 4708 _assert_(sNew>=L0); 4709 _assert_(LAB); 4710 _assert_(vertices && nbv<maxnbv); 4711 _assert_(edges && nbe<nbex); 4712 _assert_(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex); 4716 4713 4717 // new vertex on edge 4718 vb = &vertices[nbv++]; 4719 vb->m = Metric(aa,a->m,bb,b->m); 4720 vb->ReferenceNumber = e->ReferenceNumber; 4721 vb->DirOfSearch =NoDirOfSearch; 4722 double abcisse = k ? bb : aa; 4723 vb->r = e->F(abcisse); 4724 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); 4725 4726 // to take into account the direction of the edge 4727 s += lstep; 4728 edges[nbe].v[0]=va; 4729 edges[nbe].v[1]=vb; 4730 edges[nbe].ReferenceNumber =e->ReferenceNumber; 4731 edges[nbe].GeomEdgeHook = e; 4732 edges[nbe].adj[0] = PreviousNewEdge; 4733 if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe]; 4734 PreviousNewEdge=edges+nbe; 4735 nbe++; 4736 va = vb; 4714 // new vertex on edge 4715 A1=vertices+nbv++; 4716 GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge; 4717 Edge* e = edges + nbe++; 4718 double se= (sNew-L0)/LAB; 4719 if (se<0 || se>=1.000000001){ 4720 _error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]"); 4737 4721 } 4722 se = abscisseInterpole(v0.m,v1.m,AB,se,1); 4723 if (se<0 || se>1){ 4724 _error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]"); 4725 } 4726 se = k1 ? se : 1. - se; 4727 se = k1==k1equi ? se : 1. - se; 4728 VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save 4729 ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1); 4730 A1->ReferenceNumber = eeequi.ReferenceNumber; 4731 A1->DirOfSearch =NoDirOfSearch; 4732 e->GeomEdgeHook = ongequi; 4733 e->v[0]=A0; 4734 e->v[1]=A1; 4735 e->ReferenceNumber = eeequi.ReferenceNumber; 4736 e->adj[0]=PreviousNewEdge; 4738 4737 4739 /*We just added one edge to the curve: Go to the next one*/ 4740 lcurve = lcurveb; 4741 e->SetMark(); 4742 a=b; 4743 4744 /*If b is required, we are on a new curve->break*/ 4745 if (b->Required()) break; 4746 int kprev=k; 4747 k = e->AdjVertexIndex[kprev];// next vertices 4748 e = e->Adj[kprev]; 4749 _assert_(e); 4750 }// for(;;) 4751 vb = b->MeshVertexHook; 4752 4753 /*Number of edges in the last disretized curve*/ 4754 NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1); 4755 /*Number of internal vertices in the last disretized curve*/ 4756 NbNewPoints = NbEdgeCurve-1; 4757 if(!kstep){ 4758 NbVerticesOnGeomEdge0 += NbNewPoints; 4759 nbcurves++; 4738 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e; 4739 PreviousNewEdge=e; 4740 A0=A1; 4741 sNew += Lstep; 4742 if (++i== NbCreatePointOnCurve) break; 4760 4743 } 4761 nbvend=nbv+NbNewPoints;4762 lstep = lcurve / NbEdgeCurve; //approximately one4763 }// end of curve --4764 if (edges) { // last edges of the curves4765 edges[nbe].v[0]=va;4766 edges[nbe].v[1]=vb;4767 edges[nbe].ReferenceNumber = e->ReferenceNumber;4768 edges[nbe].GeomEdgeHook = e;4769 edges[nbe].adj[0] = PreviousNewEdge;4770 edges[nbe].adj[1] = 0;4771 if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe];4772 nbe++;4773 4744 } 4774 else nbe += NbEdgeCurve; 4775 } // end on curve --- 4745 4746 //some checks 4747 _assert_(ee.GeomEdgeHook->CurveNumber==ei.GeomEdgeHook->CurveNumber); 4748 if (ee[k1].GeomEdgeHook->IsRequiredVertex()) { 4749 _assert_(eeequi[k1equi].GeomEdgeHook->IsRequiredVertex()); 4750 GeomVertex * GA1 = *eeequi[k1equi].GeomEdgeHook; 4751 A1=GA1->MeshVertexHook;// the vertex in new mesh 4752 _assert_(A1-vertices>=0 && A1-vertices<nbv); 4753 break; 4754 } 4755 if (!ee.adj[k1]) { 4756 _error_("adj edge " << BTh.GetId(ee) << ", nbe=" << nbe << ", Gh.vertices=" << Gh.vertices); 4757 } 4758 pe = ee.adj[k1]; // next edge 4759 k0 = pe->Intersection(ee); 4760 peequi= eeequi.adj[k1equi]; // next edge 4761 k0equi=peequi->Intersection(eeequi); 4762 }// for(;;) end of the curve 4776 4763 } 4764 4765 if (phase){ // construction of the last edge 4766 Edge* e=edges + nbe++; 4767 e->GeomEdgeHook = ongequi; 4768 e->v[0]=A0; 4769 e->v[1]=A1; 4770 e->ReferenceNumber = peequi->ReferenceNumber; 4771 e->adj[0]=PreviousNewEdge; 4772 e->adj[1]=0; 4773 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e; 4774 PreviousNewEdge = e; 4775 4776 _assert_(i==NbCreatePointOnCurve); 4777 } 4778 4779 if (!phase) { // 4780 long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg 4781 Lstep = L/NbSegOnCurve; 4782 NbCreatePointOnCurve = NbSegOnCurve-1; 4783 NbOfNewEdge += NbSegOnCurve; 4784 NbOfNewPoints += NbCreatePointOnCurve; 4785 } 4777 4786 } 4778 } // for (i=0;i<nbe;i++) 4779 if(!step) { 4780 _assert_(!edges); 4781 _assert_(!VerticesOnGeomEdge); 4787 }// end of curve loop 4782 4788 4783 edges = new Edge[nbex=nbe]; 4784 if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0]; 4785 4786 // do the vertex on a geometrical vertex 4787 _assert_(VerticesOnGeomEdge || NbVerticesOnGeomEdge0==0); 4788 NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge; 4789 //Allocate memory 4790 if(step==0){ 4791 if(nbv+NbOfNewPoints > maxnbv) { 4792 _error_("too many vertices on geometry: " << nbv+NbOfNewPoints << " >= " << maxnbv); 4793 } 4794 edges = new Edge[NbOfNewEdge]; 4795 nbex = NbOfNewEdge; 4796 if(NbOfNewPoints) { 4797 VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints]; 4798 NbVertexOnBThEdge = NbOfNewPoints; 4799 VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints]; 4800 NbVerticesOnGeomEdgex = NbOfNewPoints; 4801 } 4802 NbOfNewPoints =0; 4803 NbOfNewEdge = 0; 4789 4804 } 4790 else{4791 _assert_(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);4792 }4793 4805 } 4806 _assert_(nbe!=0); 4807 delete [] bcurve; 4794 4808 4795 4796 4809 //Insert points inside existing triangles 4797 4810 if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n"); 4798 4811 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n");
Note:
See TracBrowser
for help on using the repository browser.
