source: issm/oecreview/Archive/18296-19100/ISSM-18492-18493.diff@ 19102

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 21.3 KB
RevLine 
[19102]1Index: ../trunk-jpl/src/c/bamg/Mesh.cpp
2===================================================================
3--- ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 18492)
4+++ ../trunk-jpl/src/c/bamg/Mesh.cpp (revision 18493)
5@@ -241,7 +241,7 @@
6 if (Gh.NbRef>0) Gh.NbRef--;
7 else if (Gh.NbRef==0) delete &Gh;
8 }
9- if(&BTh && (&BTh != this) && false){
10+ if(&BTh && (&BTh != this)){
11 if (BTh.NbRef>0) BTh.NbRef--;
12 else if (BTh.NbRef==0) delete &BTh;
13 }
14@@ -4504,295 +4504,308 @@
15 /*Get options*/
16 int verbose=bamgopts->verbose;
17
18- /*Intermediaries*/
19- int i,k;
20- int nbcurves = 0;
21- int NbNewPoints,NbEdgeCurve;
22- double lcurve,lstep,s;
23- const int MaxSubEdge = 10;
24+ Gh.NbRef++;// add a ref to Gh
25
26- R2 AB;
27- GeomVertex *a, *b;
28- BamgVertex *va,*vb;
29- GeomEdge *e;
30+ /*************************************************************************
31+ * method in 2 steps
32+ * 1 - compute the number of new edges to allocate
33+ * 2 - construct the edges
34+ * remark:
35+ * in this part we suppose to have a background mesh with the same geometry
36+ *
37+ * To construct the discretization of the new mesh we have to
38+ * rediscretize the boundary of background Mesh
39+ * because we have only the pointeur from the background mesh to the geometry.
40+ * We need the abcisse of the background mesh vertices on geometry
41+ * so a vertex is
42+ * 0 on GeomVertex ;
43+ * 1 on GeomEdge + abcisse
44+ * 2 internal
45+ *************************************************************************/
46
47- // add a ref to GH to make sure that it is not destroyed by mistake
48- Gh.NbRef++;
49+ //Check that background mesh and current mesh do have the same geometry
50+ _assert_(&BTh.Gh==&Gh);
51+ BTh.NbRef++; // add a ref to BackGround Mesh
52
53- //build background mesh flag (1 if background, else 0)
54- bool background=(&BTh != this);
55+ //Initialize new mesh
56+ BTh.SetVertexFieldOn();
57+ int* bcurve = new int[Gh.nbcurves]; //
58
59- /*Build VerticesOnGeomVertex*/
60+ /* There are 2 ways to make the loop
61+ * 1) on the geometry
62+ * 2) on the background mesh
63+ * if you do the loop on geometry, we don't have the pointeur on background,
64+ * and if you do the loop in background we have the pointeur on geometry
65+ * so do the walk on background */
66
67- //Compute the number of geometrical vertices that we are going to use to mesh
68- for (i=0;i<Gh.nbv;i++){
69- if (Gh[i].Required()) NbVerticesOnGeomVertex++;
70- }
71- //allocate
72- VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex];
73- if(NbVerticesOnGeomVertex >= maxnbv) _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
74- _assert_(nbv==0);
75- //Build VerticesOnGeomVertex
76- for (i=0;i<Gh.nbv;i++){
77- /* Add vertex only if required*/
78- if (Gh[i].Required()) {//Gh vertices Required
79+ NbVerticesOnGeomVertex=0;
80+ NbVerticesOnGeomEdge=0;
81
82- //Add the vertex
83- _assert_(nbv<maxnbv);
84- vertices[nbv]=Gh[i];
85+ /*STEP 1 copy of Required vertices*/
86
87- //Add pointer from geometry (Gh) to vertex from mesh (Th)
88- Gh[i].MeshVertexHook=vertices+nbv;
89+ int i;
90+ for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
91+ printf("\n");
92+ if(NbVerticesOnGeomVertex >= maxnbv){
93+ _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
94+ }
95
96- //Build VerticesOnGeomVertex for current point
97- VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].MeshVertexHook,Gh[i]);
98+ VerticesOnGeomVertex = new VertexOnGeom[ NbVerticesOnGeomVertex];
99+ VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex];
100
101- //nbv increment
102+ //At this point there is NO vertex but vertices should have been allocated by Init
103+ _assert_(vertices);
104+ for (i=0;i<Gh.nbv;i++){
105+ if (Gh[i].Required()) {//Gh vertices Required
106+ vertices[nbv] =Gh[i];
107+ vertices[nbv].i=I2(0,0);
108+ Gh[i].MeshVertexHook = vertices + nbv;// save Geom -> Th
109+ VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]);
110 nbv++;
111 }
112+ else Gh[i].MeshVertexHook=0;
113+ }
114+ for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){
115+ VertexOnGeom &vog=BTh.VerticesOnGeomVertex[i];
116+ if (vog.IsRequiredVertex()){
117+ GeomVertex* gv=vog;
118+ BamgVertex *bv = vog;
119+ _assert_(gv->MeshVertexHook); // use of Geom -> Th
120+ VertexOnBThVertex[NbVertexOnBThVertex++]=VertexOnVertex(gv->MeshVertexHook,bv);
121+ gv->MeshVertexHook->m = bv->m; // for taking the metric of the background mesh
122+ }
123 }
124+ _assert_(NbVertexOnBThVertex==NbVerticesOnGeomVertex); /*This might be due to MaxCornerAngle too small*/
125
126- /*Build VerticesOnGeomEdge*/
127+ /*STEP 2: reseed boundary edges*/
128
129- //check that edges is still empty (Init)
130- _assert_(!edges);
131+ // find the begining of the curve in BTh
132+ Gh.UnMarkEdges();
133+ int bfind=0;
134+ for (int i=0;i<Gh.nbcurves;i++) bcurve[i]=-1;
135
136- /* Now we are going to create the first edges corresponding
137- * to the one present in the geometry provided.
138- * We proceed in 2 steps
139- * -step 0: we count all the edges
140- * we allocate the number of edges at the end of step 0
141- * -step 1: the edges are created */
142- for (int step=0;step<2;step++){
143+ /*Loop over the backgrounf mesh BTh edges*/
144+ for (int iedge=0;iedge<BTh.nbe;iedge++){
145+ Edge &ei = BTh.edges[iedge];
146
147- //initialize number of edges and number of edges max
148- long nbex=0;
149- nbe=0;
150- long NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge;
151- Gh.UnMarkEdges();
152- nbcurves=0;
153+ /*Loop over the 2 vertices of the current edge*/
154+ for(int je=0;je<2;je++){
155
156- //go through the edges of the geometry
157- for (i=0;i<Gh.nbe;i++){
158+ /* If one of the vertex is required we are in a new curve*/
159+ if (ei[je].GeomEdgeHook->IsRequiredVertex()){
160
161- //ei = current Geometrical edge
162- GeomEdge &ei=Gh.edges[i];
163+ /*Get curve number*/
164+ int nc=ei.GeomEdgeHook->CurveNumber;
165
166- //loop over the two vertices of the edge ei
167- for(int j=0;j<2;j++) {
168+ //_printf_("Dealing with curve number " << nc << "\n");
169+ //_printf_("edge on geometry is same as GhCurve? " << (ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge)?"yes":"no\n");
170+ //if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge || ei.GeomEdgeHook==Gh.curves[nc].LastEdge){
171+ // _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");
172+ // _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");
173+ //}
174+ //BUG FIX from original bamg
175+ /*Check that we are on the same edge and right vertex (0 or 1) */
176+ if(ei.GeomEdgeHook==Gh.curves[nc].FirstEdge && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].FirstEdge)[Gh.curves[nc].FirstVertexIndex]){
177+ bcurve[nc]=iedge*2+je;
178+ bfind++;
179+ }
180+ else if ((ei.GeomEdgeHook==Gh.curves[nc].LastEdge && (GeomVertex *)*ei[je].GeomEdgeHook==&(*Gh.curves[nc].LastEdge)[Gh.curves[nc].LastVertexIndex]) && bcurve[nc]==-1){
181+ bcurve[nc]=iedge*2+je;
182+ bfind++;
183+ }
184+ }
185+ }
186+ }
187+ if (bfind!=Gh.nbcurves){
188+ delete [] bcurve;
189+ _error_("problem generating number of curves (" << Gh.nbcurves << " found in the geometry but " << bfind << " curve found in the mesh)");
190+ }
191
192- /*Take only required vertices (corner->beginning of a new curve)*/
193- if (!ei.Mark() && ei[j].Required()){
194+ // method in 2 + 1 step
195+ // 0.0) compute the length and the number of vertex to do allocation
196+ // 1.0) recompute the length
197+ // 1.1) compute the vertex
198
199- long nbvend=0;
200- Edge* PreviousNewEdge=NULL;
201- lstep = -1;
202+ long nbex=0,NbVerticesOnGeomEdgex=0;
203+ for (int step=0; step <2;step++){
204
205- /*If Edge is required (do that only once for the 2 vertices)*/
206- if(ei.Required()){
207- if (j==0){
208- //do not create internal points if required (take it as is)
209- if(step==0) nbe++;
210- else{
211- e=&ei;
212- a=ei(0);
213- b=ei(1);
214+ long NbOfNewPoints=0;
215+ long NbOfNewEdge=0;
216+ long iedge;
217+ Gh.UnMarkEdges();
218+ double L=0;
219
220- //check that edges has been allocated
221- _assert_(edges);
222- edges[nbe].v[0]=a->MeshVertexHook;
223- edges[nbe].v[1]=b->MeshVertexHook;;
224- edges[nbe].ReferenceNumber = e->ReferenceNumber;
225- edges[nbe].GeomEdgeHook = e;
226- edges[nbe].adj[0] = 0;
227- edges[nbe].adj[1] = 0;
228- nbe++;
229- }
230- }
231- }
232+ /*Go through all geometrical curve*/
233+ for (int icurve=0;icurve<Gh.nbcurves;icurve++){
234
235- /*If Edge is not required: we are on a curve*/
236- else {
237- for (int kstep=0;kstep<=step;kstep++){
238- //kstep=0, compute number of edges (discretize curve)
239- //kstep=1 create the points and edge
240- PreviousNewEdge=0;
241- NbNewPoints=0;
242- NbEdgeCurve=0;
243- if (nbvend>=maxnbv) _error_("maximum number of vertices too low! Check the domain outline or increase maxnbv");
244- lcurve =0;
245- s = lstep; //-1 initially, then length of each sub edge
246+ /*Get edge and vertex (index) of background mesh on this curve*/
247+ iedge=bcurve[icurve]/2;
248+ int jedge=bcurve[icurve]%2;
249
250- /*reminder: i = edge number, j=[0;1] vertex index in edge*/
251- k=j; // k = vertex index in edge (0 or 1)
252- e=&ei; // e = reference of current edge
253- a=ei(k); // a = pointer toward the kth vertex of the current edge
254- va = a->MeshVertexHook; // va = pointer toward mesh vertex associated
255- e->SetMark(); // Mark edge
256+ /*Get edge of Bth with index iedge*/
257+ Edge &ei = BTh.edges[iedge];
258
259- /*Loop until we reach the end of the curve*/
260- for(;;){
261- k = 1-k; // other vertx index of the curve
262- b = (*e)(k); // b = pointer toward the other vertex of the current edge
263- AB= b->r - a->r; // AB = vector of the current edge
264- Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A
265- Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B
266- double ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric
267+ /*Initialize variables*/
268+ double Lstep=0; // step between two points (phase==1)
269+ long NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1)
270
271- /* We are now creating the mesh edges from the geometrical edge selected above.
272- * The edge will be divided according to the metric previously computed and cannot
273- * be divided more than 10 times (MaxSubEdge). */
274+ /*Do phase 0 to step*/
275+ for(int phase=0;phase<=step;phase++){
276
277- //By default, there is only one subedge that is the geometrical edge itself
278- int NbSubEdge = 1;
279+ /*Current curve pointer*/
280+ Curve *curve= Gh.curves+icurve;
281
282- //initialize lSubEdge, holding the length of each subedge (cannot be higher than 10)
283- double lSubEdge[MaxSubEdge];
284+ /*Get index of current curve*/
285+ int icurveequi= Gh.GetId(curve);
286
287- //Build Subedges according to the edge length
288- if (ledge < 1.5){
289- //if ledge < 1.5 (between one and 2), take the edge as is
290- lSubEdge[0] = ledge;
291- }
292- //else, divide the edge
293- else {
294- //compute number of subedges (division of the edge), Maximum is 10
295- NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
296- /*Now, we are going to divide the edge according to the metric.
297- * Get segment by sement along the edge.
298- * Build lSubEdge, which holds the distance between the first vertex
299- * of the edge and the next point on the edge according to the
300- * discretization (each SubEdge is AB)*/
301- R2 A,B;
302- A=a->r;
303- Metric MAs=MA,MBs;
304- ledge=0;
305- double x =0, xstep= 1./NbSubEdge;
306- for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
307- x += xstep;
308- B = e->F(k ? x : 1-x);
309- MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
310- AB = A-B;
311- lSubEdge[kk]=(ledge+=(MAs(AB)+MBs(AB))/2);
312- }
313- }
314+ /*For phase 0, check that we are at the begining of the curve only*/
315+ if(phase==0 && icurveequi!=icurve) continue;
316
317- double lcurveb = lcurve+ledge;
318+ int k0=jedge,k1;
319+ Edge* pe= BTh.edges+iedge;
320+ int iedgeequi=bcurve[icurveequi]/2;
321+ int jedgeequi=bcurve[icurveequi]%2;
322
323- /*Now, create corresponding points*/
324- while(s>=lcurve && s<=lcurveb && nbv<nbvend){
325+ int k0equi=jedgeequi,k1equi;
326+ Edge * peequi= BTh.edges+iedgeequi;
327+ GeomEdge *ongequi = peequi->GeomEdgeHook;
328
329- /*Schematic of current curve
330- *
331- * a vb b // vertex
332- * 0 ll0 ll1 ledge // length from a
333- * + --- + - ... - + --S-- + --- + - ... - + // where is S
334- * 0 kk0 kk1 NbSubEdge // Sub edge index
335- *
336- */
337+ double sNew=Lstep;// abscisse of the new points (phase==1)
338+ L=0;// length of the curve
339+ long i=0;// index of new points on the curve
340+ GeomVertex * GA0 = *(*peequi)[k0equi].GeomEdgeHook;
341+ BamgVertex *A0;
342+ A0 = GA0->MeshVertexHook; // the vertex in new mesh
343+ BamgVertex *A1;
344+ VertexOnGeom *GA1;
345+ Edge* PreviousNewEdge = 0;
346
347- double ss = s-lcurve;
348+ // New Curve phase
349+ _assert_(A0-vertices>=0 && A0-vertices<nbv);
350+ if(ongequi->Required()){
351+ GeomVertex *GA1 = *(*peequi)[1-k0equi].GeomEdgeHook;
352+ A1 = GA1->MeshVertexHook; //
353+ }
354+ else {
355+ for(;;){
356+ Edge &ee=*pe;
357+ Edge &eeequi=*peequi;
358+ k1 = 1-k0; // next vertex of the edge
359+ k1equi= 1 - k0equi;
360+ _assert_(pe && ee.GeomEdgeHook);
361+ ee.GeomEdgeHook->SetMark();
362+ BamgVertex & v0=ee[0], & v1=ee[1];
363+ R2 AB=(R2)v1-(R2)v0;
364+ double L0=L,LAB;
365+ LAB=LengthInterpole(v0.m,v1.m,AB);
366+ L+= LAB;
367
368- /*Find the SubEdge containing ss using Dichotomy*/
369- int kk0=-1,kk1=NbSubEdge-1,kkk;
370- double ll0=0,ll1=ledge,llk;
371- while (kk1-kk0>1){
372- if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
373- kk1=kkk,ll1=llk;
374- else
375- kk0=kkk,ll0=llk;
376- }
377- _assert_(kk1!=kk0);
378+ if (phase){
379+ // computation of the new points for the given curve
380+ while ((i!=NbCreatePointOnCurve) && sNew<=L) {
381
382- /*Curvilinear coordinate in [0 1] of ss in current edge*/
383- // WARNING: This is what we would do
384- // ssa = (ss-ll0)/(ll1-ll0);
385- // aa = (kk0+ssa)/NbSubEdge
386- // This is what Bamg does:
387- double sbb = (ss-ll0)/(ll1-ll0);
388- /*Curvilinear coordinate in [0 1] of ss in current curve*/
389- double bb = (kk1+sbb)/NbSubEdge;
390- double aa = 1-bb;
391+ //some checks
392+ _assert_(sNew>=L0);
393+ _assert_(LAB);
394+ _assert_(vertices && nbv<maxnbv);
395+ _assert_(edges && nbe<nbex);
396+ _assert_(VerticesOnGeomEdge && NbVerticesOnGeomEdge<NbVerticesOnGeomEdgex);
397
398- // new vertex on edge
399- vb = &vertices[nbv++];
400- vb->m = Metric(aa,a->m,bb,b->m);
401- vb->ReferenceNumber = e->ReferenceNumber;
402- vb->DirOfSearch =NoDirOfSearch;
403- double abcisse = k ? bb : aa;
404- vb->r = e->F(abcisse);
405- VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);
406-
407- // to take into account the direction of the edge
408- s += lstep;
409- edges[nbe].v[0]=va;
410- edges[nbe].v[1]=vb;
411- edges[nbe].ReferenceNumber =e->ReferenceNumber;
412- edges[nbe].GeomEdgeHook = e;
413- edges[nbe].adj[0] = PreviousNewEdge;
414- if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe];
415- PreviousNewEdge=edges+nbe;
416- nbe++;
417- va = vb;
418+ // new vertex on edge
419+ A1=vertices+nbv++;
420+ GA1=VerticesOnGeomEdge+NbVerticesOnGeomEdge;
421+ Edge* e = edges + nbe++;
422+ double se= (sNew-L0)/LAB;
423+ if (se<0 || se>=1.000000001){
424+ _error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
425 }
426+ se = abscisseInterpole(v0.m,v1.m,AB,se,1);
427+ if (se<0 || se>1){
428+ _error_("Problem creating point on a boundary: se=" << se << " should be in [0 1]");
429+ }
430+ se = k1 ? se : 1. - se;
431+ se = k1==k1equi ? se : 1. - se;
432+ VertexOnBThEdge[NbVerticesOnGeomEdge++] = VertexOnEdge(A1,&eeequi,se); // save
433+ ongequi=Gh.ProjectOnCurve(eeequi,se,*A1,*GA1);
434+ A1->ReferenceNumber = eeequi.ReferenceNumber;
435+ A1->DirOfSearch =NoDirOfSearch;
436+ e->GeomEdgeHook = ongequi;
437+ e->v[0]=A0;
438+ e->v[1]=A1;
439+ e->ReferenceNumber = eeequi.ReferenceNumber;
440+ e->adj[0]=PreviousNewEdge;
441
442- /*We just added one edge to the curve: Go to the next one*/
443- lcurve = lcurveb;
444- e->SetMark();
445- a=b;
446-
447- /*If b is required, we are on a new curve->break*/
448- if (b->Required()) break;
449- int kprev=k;
450- k = e->AdjVertexIndex[kprev];// next vertices
451- e = e->Adj[kprev];
452- _assert_(e);
453- }// for(;;)
454- vb = b->MeshVertexHook;
455-
456- /*Number of edges in the last disretized curve*/
457- NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1);
458- /*Number of internal vertices in the last disretized curve*/
459- NbNewPoints = NbEdgeCurve-1;
460- if(!kstep){
461- NbVerticesOnGeomEdge0 += NbNewPoints;
462- nbcurves++;
463+ if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
464+ PreviousNewEdge=e;
465+ A0=A1;
466+ sNew += Lstep;
467+ if (++i== NbCreatePointOnCurve) break;
468 }
469- nbvend=nbv+NbNewPoints;
470- lstep = lcurve / NbEdgeCurve; //approximately one
471- }// end of curve --
472- if (edges) { // last edges of the curves
473- edges[nbe].v[0]=va;
474- edges[nbe].v[1]=vb;
475- edges[nbe].ReferenceNumber = e->ReferenceNumber;
476- edges[nbe].GeomEdgeHook = e;
477- edges[nbe].adj[0] = PreviousNewEdge;
478- edges[nbe].adj[1] = 0;
479- if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe];
480- nbe++;
481 }
482- else nbe += NbEdgeCurve;
483- } // end on curve ---
484+
485+ //some checks
486+ _assert_(ee.GeomEdgeHook->CurveNumber==ei.GeomEdgeHook->CurveNumber);
487+ if (ee[k1].GeomEdgeHook->IsRequiredVertex()) {
488+ _assert_(eeequi[k1equi].GeomEdgeHook->IsRequiredVertex());
489+ GeomVertex * GA1 = *eeequi[k1equi].GeomEdgeHook;
490+ A1=GA1->MeshVertexHook;// the vertex in new mesh
491+ _assert_(A1-vertices>=0 && A1-vertices<nbv);
492+ break;
493+ }
494+ if (!ee.adj[k1]) {
495+ _error_("adj edge " << BTh.GetId(ee) << ", nbe=" << nbe << ", Gh.vertices=" << Gh.vertices);
496+ }
497+ pe = ee.adj[k1]; // next edge
498+ k0 = pe->Intersection(ee);
499+ peequi= eeequi.adj[k1equi]; // next edge
500+ k0equi=peequi->Intersection(eeequi);
501+ }// for(;;) end of the curve
502 }
503+
504+ if (phase){ // construction of the last edge
505+ Edge* e=edges + nbe++;
506+ e->GeomEdgeHook = ongequi;
507+ e->v[0]=A0;
508+ e->v[1]=A1;
509+ e->ReferenceNumber = peequi->ReferenceNumber;
510+ e->adj[0]=PreviousNewEdge;
511+ e->adj[1]=0;
512+ if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
513+ PreviousNewEdge = e;
514+
515+ _assert_(i==NbCreatePointOnCurve);
516+ }
517+
518+ if (!phase) { //
519+ long NbSegOnCurve = Max((long)(L+0.5),(long) 1);// nb of seg
520+ Lstep = L/NbSegOnCurve;
521+ NbCreatePointOnCurve = NbSegOnCurve-1;
522+ NbOfNewEdge += NbSegOnCurve;
523+ NbOfNewPoints += NbCreatePointOnCurve;
524+ }
525 }
526- } // for (i=0;i<nbe;i++)
527- if(!step) {
528- _assert_(!edges);
529- _assert_(!VerticesOnGeomEdge);
530+ }// end of curve loop
531
532- edges = new Edge[nbex=nbe];
533- if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0];
534-
535- // do the vertex on a geometrical vertex
536- _assert_(VerticesOnGeomEdge || NbVerticesOnGeomEdge0==0);
537- NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge;
538+ //Allocate memory
539+ if(step==0){
540+ if(nbv+NbOfNewPoints > maxnbv) {
541+ _error_("too many vertices on geometry: " << nbv+NbOfNewPoints << " >= " << maxnbv);
542+ }
543+ edges = new Edge[NbOfNewEdge];
544+ nbex = NbOfNewEdge;
545+ if(NbOfNewPoints) {
546+ VerticesOnGeomEdge = new VertexOnGeom[NbOfNewPoints];
547+ NbVertexOnBThEdge = NbOfNewPoints;
548+ VertexOnBThEdge = new VertexOnEdge[NbOfNewPoints];
549+ NbVerticesOnGeomEdgex = NbOfNewPoints;
550+ }
551+ NbOfNewPoints =0;
552+ NbOfNewEdge = 0;
553 }
554- else{
555- _assert_(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0);
556- }
557 }
558+ _assert_(nbe!=0);
559+ delete [] bcurve;
560
561-
562 //Insert points inside existing triangles
563 if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n");
564 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n");
Note: See TracBrowser for help on using the repository browser.