Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/dynamic_mixed_subdivisions.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
! 2: with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
! 3: with Initial_Mixed_Cell;
! 4: with Vertices,Simplices; use Vertices,Simplices;
! 5: with Global_Dynamic_Triangulation; use Global_Dynamic_Triangulation;
! 6: with Integer_Lifting_Utilities; use Integer_Lifting_Utilities;
! 7: with Dynamic_Lifting_Functions; use Dynamic_Lifting_Functions;
! 8: with Flatten_Mixed_Subdivisions; use Flatten_Mixed_Subdivisions;
! 9: with Enumerate_Faces_of_Polytope; use Enumerate_Faces_of_Polytope;
! 10: with Common_Faces_of_Polytope; use Common_Faces_of_Polytope;
! 11: with Integer_Pruning_Methods; use Integer_Pruning_Methods;
! 12: with Mixed_Volume_Computation; use Mixed_Volume_Computation;
! 13: with Contributions_to_Mixed_Volume; use Contributions_to_Mixed_Volume;
! 14:
! 15: package body Dynamic_Mixed_Subdivisions is
! 16:
! 17: -- UTILITIES :
! 18:
! 19: function Is_Empty ( points : Array_of_Lists ) return boolean is
! 20: begin
! 21: for i in points'range loop
! 22: if not Is_Null(points(i))
! 23: then return false;
! 24: end if;
! 25: end loop;
! 26: return true;
! 27: end Is_Empty;
! 28:
! 29: function First_Non_Empty ( points : Array_of_Lists ) return integer is
! 30:
! 31: -- DESCRIPTION :
! 32: -- Returns the index of the first non empty set in the points.
! 33:
! 34: res : integer := points'first - 1;
! 35:
! 36: begin
! 37: for i in points'range loop
! 38: if not Is_Null(points(i))
! 39: then res := i;
! 40: end if;
! 41: exit when res >= points'first;
! 42: end loop;
! 43: return res;
! 44: end First_Non_Empty;
! 45:
! 46: function Is_In ( fs : Face_Structure; point : vector ) return boolean is
! 47: begin
! 48: if Is_Null(fs.t)
! 49: then return Is_In(fs.l,point);
! 50: else return Is_In(fs.t,point);
! 51: end if;
! 52: end Is_In;
! 53:
! 54: -- procedure put ( n : in natural; fs : in Face_Structure ) is
! 55: -- begin
! 56: -- put_line("The list of lifted points : "); put(fs.l);
! 57: -- put_line("The list of k-faces : "); put(fs.f);
! 58: -- put_line("The triangulation : "); put(n,fs.t);
! 59: -- end put;
! 60:
! 61: -- procedure put ( n : in natural; fs : in Face_Structures ) is
! 62: -- begin
! 63: -- for i in fs'range loop
! 64: -- put("face structure for component "); put(i,1);
! 65: -- put_line(" :");
! 66: -- put(n,fs(i));
! 67: -- end loop;
! 68: -- end put;
! 69:
! 70: -- FLATTENING :
! 71:
! 72: procedure Flatten ( points : in out VecVec ) is
! 73:
! 74: pt : Link_to_Vector;
! 75:
! 76: begin
! 77: for i in points'range loop
! 78: pt := points(i);
! 79: if pt(pt'last) /= 0
! 80: then pt(pt'last) := 0;
! 81: end if;
! 82: end loop;
! 83: end Flatten;
! 84:
! 85: procedure Flatten ( f : in out Face ) is
! 86: begin
! 87: Flatten(f.all);
! 88: end Flatten;
! 89:
! 90: procedure Flatten ( fs : in out Faces ) is
! 91:
! 92: tmp : Faces := fs;
! 93: f : Face;
! 94:
! 95: begin
! 96: while not Is_Null(tmp) loop
! 97: f := Head_Of(tmp);
! 98: Flatten(f);
! 99: Set_Head(tmp,f);
! 100: tmp := Tail_Of(tmp);
! 101: end loop;
! 102: end Flatten;
! 103:
! 104: procedure Flatten ( fs : in out Face_Structure ) is
! 105: begin
! 106: Flatten(fs.l); Flatten(fs.t); Flatten(fs.f);
! 107: end Flatten;
! 108:
! 109: procedure Flatten ( fs : in out Face_Structures ) is
! 110: begin
! 111: for i in fs'range loop
! 112: Flatten(fs(i));
! 113: end loop;
! 114: end Flatten;
! 115:
! 116: -- ZERO CONTRIBUTION :
! 117:
! 118: function Extract_Facet ( n : natural; s : Simplex ) return Face is
! 119:
! 120: ver : constant VecVec := Simplices.Vertices(s);
! 121: res : Face := new VecVec(ver'range);
! 122:
! 123: begin
! 124: for i in res'range loop
! 125: res(i) := new Standard_Integer_Vectors.Vector'(ver(i)(1..n));
! 126: end loop;
! 127: return res;
! 128: end Extract_Facet;
! 129:
! 130: function Extract_Facets ( n : natural; t : Triangulation; x : Vector )
! 131: return Faces is
! 132:
! 133: -- DESCRIPTION :
! 134: -- Returns the list of facets in t that all contain x.
! 135: -- The facets are given in their original coordinates.
! 136:
! 137: -- REQUIRED : x is lifted, i.e., x'length = n+1.
! 138:
! 139: tmp : Triangulation := t;
! 140: res,res_last : Faces;
! 141:
! 142: begin
! 143: while not Is_Null(tmp) loop
! 144: declare
! 145: s : Simplex := Head_Of(tmp);
! 146: begin
! 147: if Is_Vertex(s,x)
! 148: then Append(res,res_last,Extract_Facet(n,s));
! 149: end if;
! 150: end;
! 151: tmp := Tail_Of(tmp);
! 152: end loop;
! 153: return res;
! 154: end Extract_Facets;
! 155:
! 156: function Zero_Contribution ( n : natural; fs : Face_Structures;
! 157: x : Vector; i : natural ) return boolean is
! 158:
! 159: -- DESCRIPTION :
! 160: -- Returns true if the point x does not contribute to the mixed volume
! 161: -- when considered to the ith component of the already lifted points.
! 162:
! 163: -- REQUIRED : x'length = n+1, n = dimension before lifting.
! 164:
! 165: res : boolean := false;
! 166: supp : Array_of_Lists(fs'range);
! 167:
! 168: begin
! 169: for j in supp'range loop
! 170: supp(j) := Reduce(fs(j).l,n+1);
! 171: end loop;
! 172: if Length_Of(fs(i).t) > 1
! 173: then declare
! 174: facets : Faces := Extract_Facets(n,fs(i).t,x);
! 175: begin
! 176: res := Simple_Zero_Contribution(supp,facets,x(1..n),i);
! 177: end;
! 178: else res := Simple_Zero_Contribution(supp,x(1..n),i);
! 179: end if;
! 180: Deep_Clear(supp);
! 181: return res;
! 182: end Zero_Contribution;
! 183:
! 184: -- INITIALIZATION :
! 185:
! 186: procedure Initialize
! 187: ( n : in natural; mix : in Vector; points : in Array_of_Lists;
! 188: mixsub : in out Mixed_Subdivision;
! 189: fs : in out Face_Structures; rest : in out Array_of_Lists;
! 190: done : in out boolean ) is
! 191:
! 192: -- DESCRIPTION :
! 193: -- Performs initialization of the dynamic lifting algorithm.
! 194:
! 195: -- ON ENTRY :
! 196: -- n length of the vectors in points;
! 197: -- mix type of mixture;
! 198: -- mixsub mixed subdivision of the lifted points;
! 199: -- fs face structures of the lifted points.
! 200:
! 201: -- ON RETURN :
! 202: -- mixsub if empty on entry, then it contains the initial mixed
! 203: -- cell, in case the problem is nondegenerate;
! 204: -- rest rest of point list to be processed;
! 205: -- done if true, then either the mixed volume equals zero,
! 206: -- or there are no more points to process.
! 207:
! 208: null_points : Array_of_Lists(points'range);
! 209:
! 210: begin
! 211: if Is_Null(mixsub)
! 212: then -- compute an initial mixed cell
! 213: declare
! 214: mic : Mixed_Cell;
! 215: begin
! 216: Initial_Mixed_Cell(n,mix,points,mic,rest);
! 217: if Mixed_Volume(n,mix,mic) > 0 -- check on degeneracy
! 218: then
! 219: -- initialize the mixed subdivision :
! 220: Construct(mic,mixsub);
! 221: -- initialize the face structures :
! 222: for i in mic.pts'range loop
! 223: declare
! 224: tmp : List := mic.pts(i);
! 225: fc : Face;
! 226: begin
! 227: while not Is_Null(tmp) loop
! 228: Append(fs(i).l,fs(i).last,Head_Of(tmp).all);
! 229: tmp := Tail_of(tmp);
! 230: end loop;
! 231: fc := new VecVec'(Shallow_Create(mic.pts(i)));
! 232: Construct(fc,fs(i).f);
! 233: end;
! 234: end loop;
! 235: if mix'last = mix'first
! 236: then declare
! 237: v : constant VecVec := Shallow_Create(fs(fs'first).l);
! 238: s : Simplex := Create(v);
! 239: begin
! 240: fs(fs'first).t := Create(s);
! 241: end;
! 242: end if;
! 243: done := Is_Empty(rest);
! 244: else -- degenerate problem => mixed volume = 0
! 245: rest := null_points; done := true;
! 246: end if;
! 247: end;
! 248: else
! 249: rest := points;
! 250: done := Is_Empty(rest);
! 251: end if;
! 252: end Initialize;
! 253:
! 254: function Initial_Triangulation
! 255: ( n : in natural; l : in List; point : in Link_to_Vector )
! 256: return Triangulation is
! 257:
! 258: -- DESCRIPTION :
! 259: -- Given a list of lifted points with Length_Of(l) >= n and another
! 260: -- lifted point, a nonempty triangulation will be returned, when the
! 261: -- volume of \conv(l U {point}) > 0.
! 262:
! 263: res : Triangulation;
! 264: del,span : List;
! 265: delpoint : Link_to_Vector := new vector(1..n);
! 266:
! 267: function Deflate ( lifted : in List ) return List is
! 268:
! 269: -- DESCRIPTION :
! 270: -- Discards the lifting value of the points in the list l.
! 271:
! 272: res,res_last : List;
! 273: tmp : List := lifted;
! 274: pt : Link_to_Vector;
! 275:
! 276: begin
! 277: while not Is_Null(tmp) loop
! 278: pt := Head_Of(tmp);
! 279: declare
! 280: dpt : Link_to_Vector := new vector(1..n);
! 281: begin
! 282: dpt.all := pt(dpt'range);
! 283: Append(res,res_last,dpt);
! 284: end;
! 285: tmp := Tail_Of(tmp);
! 286: end loop;
! 287: return res;
! 288: end Deflate;
! 289:
! 290: function Lift ( pt,liftpt : Link_to_Vector; lifted : List )
! 291: return Link_to_Vector is
! 292:
! 293: -- DESCRIPTION :
! 294: -- Compares the value of pt with that of liftpt and
! 295: -- looks up the lifting value of the point pt in the list lifted.
! 296: -- The lifted point will be returned.
! 297:
! 298: -- REQUIRED :
! 299: -- Either the point pt should be equal to Deflate(liftpt)
! 300: -- or it should belong to Deflate(lifted).
! 301:
! 302: tmp : List := lifted;
! 303: res : Link_to_Vector;
! 304:
! 305: begin
! 306: if pt.all = liftpt(pt'range)
! 307: then return liftpt;
! 308: else while not Is_Null(tmp) loop
! 309: res := Head_Of(tmp);
! 310: if pt.all = res(pt'range)
! 311: then return res;
! 312: else tmp := Tail_Of(tmp);
! 313: end if;
! 314: end loop;
! 315: return res;
! 316: end if;
! 317: end Lift;
! 318:
! 319: begin
! 320: del := Deflate(l);
! 321: delpoint.all := point(delpoint'range);
! 322: Construct(delpoint,del);
! 323: span := Extremal_Points(n,del);
! 324: if Length_Of(span) = n+1
! 325: then
! 326: declare
! 327: verpts : VecVec(1..n+1);
! 328: tmp : List := span;
! 329: pt : Link_to_Vector;
! 330: s : Simplex;
! 331: begin
! 332: for i in verpts'range loop
! 333: verpts(i) := Lift(Head_Of(tmp),point,l);
! 334: tmp := Tail_Of(tmp);
! 335: end loop;
! 336: s := Create(verpts);
! 337: res := Create(s);
! 338: tmp := l;
! 339: while not Is_Null(tmp) loop
! 340: pt := Head_Of(tmp);
! 341: if not Is_In(span,pt)
! 342: then Update(res,pt);
! 343: end if;
! 344: tmp := Tail_Of(tmp);
! 345: end loop;
! 346: end;
! 347: end if;
! 348: Clear(del); Clear(span); Clear(delpoint);
! 349: return res;
! 350: end Initial_Triangulation;
! 351:
! 352: -- CHOOSE NEXT POINT :
! 353:
! 354: procedure Next_Point
! 355: ( n : in natural; points : in out Array_of_Lists;
! 356: order : in boolean;
! 357: index : in out integer; point : out Link_to_Vector ) is
! 358:
! 359: -- DESCRIPTION :
! 360: -- Chooses a point out of the lists of points.
! 361:
! 362: -- ON ENTRY :
! 363: -- n length of the elements in the point lists;
! 364: -- points nonempty lists of points;
! 365: -- order if true, then the first element out of the next first
! 366: -- nonempty list will be chosen;
! 367: -- index indicates component where not Is_Null(points(index)).
! 368:
! 369: -- ON RETURN :
! 370: -- points lists of points with the point removed;
! 371: -- index points to the next nonempty list,
! 372: -- if index < points'first, then Is_Empty(points);
! 373: -- point the next point to be processed.
! 374:
! 375: newindex : integer := points'first-1;
! 376: pt : Link_to_Vector;
! 377:
! 378: begin
! 379: if order
! 380: then pt := Head_Of(points(index));
! 381: else pt := Max_Extreme(points(index),n,-3,3);
! 382: Swap_to_Front(points(index),pt);
! 383: end if;
! 384: points(index) := Tail_Of(points(index));
! 385: for i in (index+1)..points'last loop
! 386: if not Is_Null(points(i))
! 387: then newindex := i;
! 388: end if;
! 389: exit when newindex >= points'first;
! 390: end loop;
! 391: if newindex < points'first
! 392: then for i in points'first..index loop
! 393: if not Is_Null(points(i))
! 394: then newindex := i;
! 395: end if;
! 396: exit when newindex >= points'first;
! 397: end loop;
! 398: end if;
! 399: index := newindex;
! 400: point := pt;
! 401: end Next_Point;
! 402:
! 403: -- UPDATE ROUTINES :
! 404:
! 405: procedure Merge ( mic : in out Mixed_Cell;
! 406: mixsub : in out Mixed_Subdivision ) is
! 407:
! 408: tmp : Mixed_Subdivision := mixsub;
! 409: done : boolean := false;
! 410:
! 411: begin
! 412: while not Is_Null(tmp) loop
! 413: declare
! 414: mic1 : Mixed_Cell := Head_Of(tmp);
! 415: last : List;
! 416: begin
! 417: if Equal(mic.nor,mic1.nor)
! 418: then for k in mic1.pts'range loop
! 419: last := mic1.pts(k);
! 420: while not Is_Null(Tail_Of(last)) loop
! 421: last := Tail_Of(last);
! 422: end loop;
! 423: Deep_Concat_Diff(mic.pts(k),mic1.pts(k),last);
! 424: end loop;
! 425: done := true;
! 426: else tmp := Tail_Of(tmp);
! 427: end if;
! 428: end;
! 429: exit when done;
! 430: end loop;
! 431: if done
! 432: then Deep_Clear(mic);
! 433: else Construct(mic,mixsub);
! 434: end if;
! 435: end Merge;
! 436:
! 437: procedure Merge ( cells : in Mixed_Subdivision;
! 438: mixsub : in out Mixed_Subdivision ) is
! 439:
! 440: tmp : Mixed_Subdivision := cells;
! 441:
! 442: begin
! 443: while not Is_Null(tmp) loop
! 444: declare
! 445: mic : Mixed_Cell := Head_Of(tmp);
! 446: begin
! 447: Merge(mic,mixsub);
! 448: end;
! 449: tmp := Tail_Of(tmp);
! 450: end loop;
! 451: end Merge;
! 452:
! 453: procedure Compute_New_Faces
! 454: ( fs : in out Face_Structure; k,n : in natural;
! 455: point : in out Link_to_Vector; newfs : out Faces ) is
! 456:
! 457: -- DESCRIPTION :
! 458: -- Given a point, the new faces will be computed and returned.
! 459:
! 460: -- ON ENTRY :
! 461: -- fs a face structure;
! 462: -- k dimension of the faces to generate;
! 463: -- n dimension without the lifting;
! 464: -- point point that has to be in all the faces.
! 465:
! 466: -- ON RETURN :
! 467: -- fs face structure with updated triangulation,
! 468: -- the new faces are not added, also the list is not updated;
! 469: -- point lifted conservatively w.r.t. fs.t;
! 470: -- newfs k-faces, which all contain the given point.
! 471:
! 472: res,res_last : Faces;
! 473:
! 474: procedure Append ( fc : in VecVec ) is
! 475:
! 476: f : Face;
! 477:
! 478: begin
! 479: f := new VecVec'(fc);
! 480: Append(res,res_last,f);
! 481: end Append;
! 482: procedure EnumLis is new Enumerate_Lower_Faces_in_List(Append);
! 483: procedure EnumTri is new Enumerate_Faces_in_Triangulation(Append);
! 484:
! 485: begin
! 486: -- COMPUTE THE NEW FACES AND UPDATE fs :
! 487: if Is_Null(fs.t)
! 488: then
! 489: if Length_Of(fs.l) >= n
! 490: then
! 491: fs.t := Initial_Triangulation(n,fs.l,point);
! 492: if Is_Null(fs.t)
! 493: then EnumLis(fs.l,point,k);
! 494: else EnumTri(fs.t,point,k);
! 495: end if;
! 496: else
! 497: EnumLis(fs.l,point,k);
! 498: end if;
! 499: else
! 500: declare
! 501: newt : Triangulation;
! 502: begin
! 503: point(point'last) := Lift_to_Place(fs.t,point.all);
! 504: Update(fs.t,point,newt);
! 505: Enumtri(newt,point,k);
! 506: end;
! 507: end if;
! 508: Append(fs.l,fs.last,point);
! 509: newfs := res;
! 510: end Compute_New_Faces;
! 511:
! 512: procedure Swap ( index : in natural; points : in out Array_of_Lists ) is
! 513:
! 514: -- DESCRIPTION :
! 515: -- Swaps the first component of points with the component index.
! 516:
! 517: tmp : List := points(index);
! 518:
! 519: begin
! 520: points(index) := points(points'first);
! 521: points(points'first) := tmp;
! 522: end Swap;
! 523:
! 524: procedure Swap ( index : in natural; mixsub : in out Mixed_Subdivision ) is
! 525:
! 526: -- DESCRIPTION :
! 527: -- Swaps the first component of each cell with the component index.
! 528:
! 529: tmp : Mixed_Subdivision := mixsub;
! 530:
! 531: begin
! 532: while not Is_Null(tmp) loop
! 533: declare
! 534: mic : Mixed_Cell := Head_Of(tmp);
! 535: begin
! 536: Swap(index,mic.pts.all);
! 537: -- Set_Head(tmp,mic);
! 538: end;
! 539: tmp := Tail_Of(tmp);
! 540: end loop;
! 541: end Swap;
! 542:
! 543: procedure Create_Cells
! 544: ( index,n : in natural; mix : in Vector;
! 545: afs : in Array_of_Faces; lif : in Array_of_Lists;
! 546: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
! 547: mixsub : out Mixed_Subdivision ) is
! 548:
! 549: -- DESCRIPTION :
! 550: -- Creates a mixed subdivision by considering the faces
! 551: -- of component index first.
! 552:
! 553: res : Mixed_Subdivision;
! 554: wrkmix : Vector(mix'range) := mix;
! 555: wrkafs : Array_of_Faces(afs'range) := afs;
! 556: wrklif : Array_of_Lists(lif'range) := lif;
! 557:
! 558: begin
! 559: wrkmix(wrkmix'first) := mix(index); wrkmix(index) := mix(mix'first);
! 560: wrkafs(wrkafs'first) := afs(index); wrkafs(index) := afs(afs'first);
! 561: wrklif(wrklif'first) := lif(index); wrklif(index) := lif(lif'first);
! 562: Create_CS(n,wrkmix,wrkafs,wrklif,nbsucc,nbfail,res);
! 563: Swap(index,res);
! 564: mixsub := res;
! 565: end Create_Cells;
! 566:
! 567: procedure Compute_New_Cells
! 568: ( n : in natural; mix : in Vector; mic : in Mixed_Cell;
! 569: afs : in Array_of_Faces; index : in integer;
! 570: lifted : in Array_of_Lists; mixsub : out Mixed_Subdivision;
! 571: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
! 572:
! 573: -- DESCRIPTION :
! 574: -- Computes the new cells by considering a given mixed cell with
! 575: -- a tuple of faces.
! 576:
! 577: -- ON ENTRY :
! 578: -- n dimension before lifting;
! 579: -- mix type of mixture;
! 580: -- mic a given mixed cell;
! 581: -- afs type of faces;
! 582: -- index component where new faces are computed for;
! 583: -- lifted tuple of lifted points;
! 584: -- nbsucc number of succesful tests of face-face combinations;
! 585: -- nbfail number of unsuccesful tests of face-face combinations.
! 586:
! 587: -- ON RETURN :
! 588: -- mixsub the new mixed cells;
! 589: -- nbsucc updated number of succesful tests;
! 590: -- nbfail updated number of unsuccesful tests;
! 591:
! 592: res,res2 : Mixed_Subdivision;
! 593: neiafs : Array_of_Faces(afs'range);
! 594:
! 595: begin
! 596: neiafs(index) := Neighboring_Faces(mic,afs(index),index);
! 597: if not Is_Null(neiafs(index))
! 598: then
! 599: for i in neiafs'range loop
! 600: if i /= index
! 601: then neiafs(i) := Neighboring_Faces(mic,afs(i),i);
! 602: end if;
! 603: end loop;
! 604: if index /= 1
! 605: then Create_Cells(index,n,mix,neiafs,lifted,nbsucc,nbfail,res);
! 606: else Create_CS(n,mix,neiafs,lifted,nbsucc,nbfail,res);
! 607: end if;
! 608: for i in neiafs'range loop
! 609: if i /= index
! 610: then Shallow_Clear(neiafs(i));
! 611: end if;
! 612: end loop;
! 613: end if;
! 614: Shallow_Clear(neiafs(index));
! 615: res2 := Create(lifted,res); -- test on normals
! 616: Shallow_Clear(res);
! 617: mixsub := res2;
! 618: end Compute_New_Cells;
! 619:
! 620: procedure Compute_New_Cells
! 621: ( n : in natural; mix : in Vector;
! 622: mixsub : in Mixed_Subdivision; index : in integer;
! 623: newfaces : in Faces; fs : in Face_Structures;
! 624: newcells : out Mixed_Subdivision;
! 625: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
! 626:
! 627: res : Mixed_Subdivision;
! 628: afs : Array_of_Faces(fs'range);
! 629: lifted : Array_of_Lists(fs'range);
! 630:
! 631: begin
! 632: -- CREATE NEW CELLS ONLY WITH THE NEW FACES :
! 633: for i in afs'range loop
! 634: if i = index
! 635: then afs(index) := newfaces;
! 636: else afs(i) := fs(i).f;
! 637: end if;
! 638: lifted(i) := fs(i).l;
! 639: end loop;
! 640: if not Is_Null(fs(index).t)
! 641: then
! 642: -- ENUMERATE ALL CELLS IN mixsub :
! 643: declare
! 644: tmp : Mixed_Subdivision := mixsub;
! 645: begin
! 646: while not Is_Null(tmp) loop
! 647: declare
! 648: newcells2 : Mixed_Subdivision;
! 649: begin
! 650: Compute_New_Cells(n,mix,Head_Of(tmp),afs,index,lifted,
! 651: newcells2,nbsucc,nbfail);
! 652: Merge(newcells2,res); Shallow_Clear(newcells2);
! 653: end;
! 654: tmp := Tail_Of(tmp);
! 655: end loop;
! 656: end;
! 657: else
! 658: -- COMPUTE ALL NEW CELLS AT ONCE :
! 659: declare
! 660: res2 : Mixed_Subdivision;
! 661: begin
! 662: if index /= 1
! 663: then Create_Cells(index,n,mix,afs,lifted,nbsucc,nbfail,res2);
! 664: else Create_CS(n,mix,afs,lifted,nbsucc,nbfail,res2);
! 665: end if;
! 666: res := Create(lifted,res2);
! 667: Shallow_Clear(res2);
! 668: end;
! 669: end if;
! 670: newcells := res;
! 671: end Compute_New_Cells;
! 672:
! 673: -- BASIC VERSION : WITHOUT OUTPUT GENERICS :
! 674:
! 675: procedure Dynamic_Lifting
! 676: ( n : in natural; mix : in Vector;
! 677: points : in Array_of_Lists;
! 678: order,inter,conmv : in boolean; maxli : in natural;
! 679: mixsub : in out Mixed_Subdivision;
! 680: fs : in out Face_Structures;
! 681: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
! 682:
! 683: rest,inner : Array_of_Lists(points'range);
! 684: index,newindex : integer;
! 685: finished : boolean := false;
! 686: pt : Link_to_Vector;
! 687: nexli : natural := 1;
! 688:
! 689: begin
! 690: Initialize(n,mix,points,mixsub,fs,rest,finished);
! 691: if not finished
! 692: then
! 693: index := First_Non_Empty(rest); newindex := index;
! 694: while not finished loop -- ITERATE UNTIL NO MORE POINTS :
! 695: Next_Point(n,rest,order,newindex,pt);
! 696: declare
! 697: point : Link_to_Vector := new vector(pt'first..pt'last+1);
! 698: newfa : Faces;
! 699: newcells : Mixed_Subdivision;
! 700: begin -- LIFT THE POINT CONSERVATIVELY :
! 701: point(pt'range) := pt.all;
! 702: point(point'last) := 1;
! 703: if inter and then Is_In(fs(index),point.all)
! 704: then
! 705: Clear(point); Construct(pt,inner(index));
! 706: else
! 707: nexli := Conservative_Lifting(mixsub,index,point.all);
! 708: if (maxli > 0) and then nexli > maxli
! 709: then Flatten(mixsub); Flatten(fs);
! 710: nexli := 1;
! 711: end if;
! 712: point(point'last) := nexli;
! 713: -- COMPUTE NEW FACES AND NEW CELLS :
! 714: Compute_New_Faces(fs(index),mix(index),n,point,newfa);
! 715: if not conmv or else not Zero_Contribution(n,fs,point.all,index)
! 716: then Compute_New_Cells(n,mix,mixsub,index,newfa,fs,
! 717: newcells,nbsucc,nbfail);
! 718: -- UPDATE THE MIXED SUBDIVISION AND THE FACE STRUCTURES :
! 719: Construct(newcells,mixsub); Shallow_Clear(newcells);
! 720: end if;
! 721: Construct(fs(index).f,newfa); Shallow_Clear(newfa);
! 722: end if;
! 723: end;
! 724: index := newindex;
! 725: finished := (index < rest'first);
! 726: end loop;
! 727: end if;
! 728: if inter -- lift out the interior points
! 729: then for i in inner'range loop
! 730: Constant_Lifting(inner(i),nexli,fs(i).l,fs(i).last);
! 731: Shallow_Clear(inner(i));
! 732: end loop;
! 733: end if;
! 734: exception
! 735: when numeric_error -- probably due to a too high lifting
! 736: => Flatten(mixsub); Flatten(fs);
! 737: Dynamic_Lifting(n,mix,rest,order,inter,conmv,maxli,mixsub,fs,
! 738: nbsucc,nbfail);
! 739: end Dynamic_Lifting;
! 740:
! 741: -- EXTENDED VERSIONS : WITH OUTPUT GENERICS :
! 742:
! 743: procedure Dynamic_Lifting_with_Flat
! 744: ( n : in natural; mix : in Vector; points : in Array_of_Lists;
! 745: order,inter,conmv : in boolean; maxli : in natural;
! 746: mixsub : in out Mixed_Subdivision;
! 747: fs : in out Face_Structures;
! 748: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
! 749:
! 750: rest,inner : Array_of_Lists(points'range);
! 751: index,newindex : integer;
! 752: finished : boolean := false;
! 753: pt : Link_to_Vector;
! 754: nexli : natural := 1;
! 755:
! 756: begin
! 757: Initialize(n,mix,points,mixsub,fs,rest,finished);
! 758: if not finished
! 759: then
! 760: index := First_Non_Empty(rest); newindex := index;
! 761: while not finished loop -- ITERATE UNTIL NO MORE POINTS :
! 762: Next_Point(n,rest,order,newindex,pt);
! 763: declare
! 764: point : Link_to_Vector := new vector(pt'first..pt'last+1);
! 765: newfa : Faces;
! 766: newcells : Mixed_Subdivision;
! 767: begin -- LIFT THE POINT CONSERVATIVELY :
! 768: point(pt'range) := pt.all;
! 769: point(point'last) := 1;
! 770: if inter and then Is_In(fs(index),point.all)
! 771: then
! 772: Clear(point); Construct(pt,inner(index));
! 773: else
! 774: nexli := Conservative_Lifting(mixsub,index,point.all);
! 775: if (maxli > 0) and then nexli > maxli
! 776: then Before_Flattening(mixsub,fs);
! 777: Flatten(mixsub); Flatten(fs);
! 778: nexli := 1;
! 779: end if;
! 780: point(point'last) := nexli;
! 781: -- COMPUTE NEW FACES AND NEW CELLS :
! 782: Compute_New_Faces(fs(index),mix(index),n,point,newfa);
! 783: if not conmv or else not Zero_Contribution(n,fs,point.all,index)
! 784: then Compute_New_Cells(n,mix,mixsub,index,newfa,fs,
! 785: newcells,nbsucc,nbfail);
! 786: -- UPDATE THE MIXED SUBDIVISION AND THE FACE STRUCTURES :
! 787: Construct(newcells,mixsub); Shallow_Clear(newcells);
! 788: end if;
! 789: Construct(fs(index).f,newfa); Shallow_Clear(newfa);
! 790: end if;
! 791: end;
! 792: index := newindex;
! 793: finished := (index < rest'first);
! 794: end loop;
! 795: end if;
! 796: if inter -- lift out the interior points
! 797: then for i in inner'range loop
! 798: Constant_Lifting(inner(i),nexli,fs(i).l,fs(i).last);
! 799: Shallow_Clear(inner(i));
! 800: end loop;
! 801: end if;
! 802: exception
! 803: when numeric_error -- probably due to a too high lifting
! 804: => Before_Flattening(mixsub,fs);
! 805: Flatten(mixsub); Flatten(fs);
! 806: Dynamic_Lifting_with_Flat(n,mix,rest,order,inter,conmv,maxli,mixsub,
! 807: fs,nbsucc,nbfail);
! 808: end Dynamic_Lifting_with_Flat;
! 809:
! 810: procedure Dynamic_Lifting_with_New
! 811: ( n : in natural; mix : in Vector; points : in Array_of_Lists;
! 812: order,inter,conmv : in boolean; maxli : in natural;
! 813: mixsub : in out Mixed_Subdivision;
! 814: fs : in out Face_Structures;
! 815: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
! 816:
! 817: rest,inner : Array_of_Lists(points'range);
! 818: index,newindex : integer;
! 819: finished : boolean := false;
! 820: pt : Link_to_Vector;
! 821: nexli : natural := 1;
! 822:
! 823: begin
! 824: Initialize(n,mix,points,mixsub,fs,rest,finished);
! 825: Process_New_Cells(mixsub,0,Head_Of(fs(fs'last).last).all);
! 826: if not finished
! 827: then
! 828: index := First_Non_Empty(rest); newindex := index;
! 829: while not finished loop -- ITERATE UNTIL NO MORE POINTS :
! 830: Next_Point(n,rest,order,newindex,pt);
! 831: declare
! 832: point : Link_to_Vector := new vector(pt'first..pt'last+1);
! 833: newfa : Faces;
! 834: newcells : Mixed_Subdivision;
! 835: begin -- LIFT THE POINT CONSERVATIVELY :
! 836: point(pt'range) := pt.all;
! 837: point(point'last) := 1;
! 838: if inter and then Is_In(fs(index),point.all)
! 839: then
! 840: Clear(point); Construct(pt,inner(index));
! 841: else
! 842: nexli := Conservative_Lifting(mixsub,index,point.all);
! 843: if (maxli > 0) and then nexli > maxli
! 844: then Flatten(mixsub); Flatten(fs);
! 845: nexli := 1;
! 846: end if;
! 847: point(point'last) := nexli;
! 848: -- COMPUTE NEW FACES AND NEW CELLS :
! 849: Compute_New_Faces(fs(index),mix(index),n,point,newfa);
! 850: if not conmv or else not Zero_Contribution(n,fs,point.all,index)
! 851: then Compute_New_Cells(n,mix,mixsub,index,newfa,fs,newcells,
! 852: nbsucc,nbfail);
! 853: Process_New_Cells(newcells,index,point.all);
! 854: -- UPDATE THE MIXED SUBDIVISION AND THE FACE STRUCTURES :
! 855: Construct(newcells,mixsub); Shallow_Clear(newcells);
! 856: end if;
! 857: Construct(fs(index).f,newfa); Shallow_Clear(newfa);
! 858: end if;
! 859: end;
! 860: index := newindex;
! 861: finished := (index < rest'first);
! 862: end loop;
! 863: end if;
! 864: if inter -- lift out the interior points
! 865: then for i in inner'range loop
! 866: Constant_Lifting(inner(i),nexli,fs(i).l,fs(i).last);
! 867: Shallow_Clear(inner(i));
! 868: end loop;
! 869: end if;
! 870: exception
! 871: when numeric_error -- probably due to a too high lifting
! 872: => Flatten(mixsub); Flatten(fs);
! 873: Dynamic_Lifting_with_New(n,mix,rest,order,inter,conmv,maxli,mixsub,fs,
! 874: nbsucc,nbfail);
! 875: end Dynamic_Lifting_with_New;
! 876:
! 877: procedure Dynamic_Lifting_with_Flat_and_New
! 878: ( n : in natural; mix : in Vector; points : in Array_of_Lists;
! 879: order,inter,conmv : in boolean; maxli : in natural;
! 880: mixsub : in out Mixed_Subdivision;
! 881: fs : in out Face_Structures;
! 882: nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
! 883:
! 884: rest,inner : Array_of_Lists(points'range);
! 885: index,newindex : integer;
! 886: finished : boolean := false;
! 887: pt : Link_to_Vector;
! 888: nexli : natural := 1;
! 889:
! 890: begin
! 891: Initialize(n,mix,points,mixsub,fs,rest,finished);
! 892: Process_New_Cells(mixsub,0,Head_Of(fs(fs'last).last).all);
! 893: if not finished
! 894: then
! 895: index := First_Non_Empty(rest); newindex := index;
! 896: while not finished loop -- ITERATE UNTIL NO MORE POINTS
! 897: Next_Point(n,rest,order,newindex,pt);
! 898: declare
! 899: point : Link_to_Vector := new vector(pt'first..pt'last+1);
! 900: newfa : Faces;
! 901: newcells : Mixed_Subdivision;
! 902: begin -- LIFT THE POINT CONSERVATIVELY :
! 903: point(pt'range) := pt.all;
! 904: point(point'last) := 1;
! 905: if inter and then Is_In(fs(index),point.all)
! 906: then
! 907: Clear(point); Construct(pt,inner(index));
! 908: else
! 909: nexli := Conservative_Lifting(mixsub,index,point.all);
! 910: if (maxli > 0) and then nexli > maxli
! 911: then Before_Flattening(mixsub,fs);
! 912: Flatten(mixsub); Flatten(fs);
! 913: nexli := 1;
! 914: end if;
! 915: point(point'last) := nexli;
! 916: -- COMPUTE NEW FACES AND NEW CELLS :
! 917: Compute_New_Faces(fs(index),mix(index),n,point,newfa);
! 918: if not conmv or else not Zero_Contribution(n,fs,point.all,index)
! 919: then Compute_New_Cells(n,mix,mixsub,index,newfa,fs,newcells,
! 920: nbsucc,nbfail);
! 921: Process_New_Cells(newcells,index,point.all);
! 922: -- UPDATE THE MIXED SUBDIVISION AND THE FACE STRUCTURES :
! 923: Construct(newcells,mixsub); Shallow_Clear(newcells);
! 924: end if;
! 925: Construct(fs(index).f,newfa); Shallow_Clear(newfa);
! 926: end if;
! 927: end;
! 928: index := newindex;
! 929: finished := (index < rest'first);
! 930: end loop;
! 931: end if;
! 932: if inter -- lift out the interior points
! 933: then for i in inner'range loop
! 934: Constant_Lifting(inner(i),nexli,fs(i).l,fs(i).last);
! 935: Shallow_Clear(inner(i));
! 936: end loop;
! 937: end if;
! 938: exception
! 939: when numeric_error -- probably due to a too high lifting
! 940: => Before_Flattening(mixsub,fs);
! 941: Flatten(mixsub); Flatten(fs);
! 942: Dynamic_Lifting_with_Flat_and_New
! 943: (n,mix,rest,order,inter,conmv,maxli,mixsub,fs,nbsucc,nbfail);
! 944: end Dynamic_Lifting_with_Flat_and_New;
! 945:
! 946: end Dynamic_Mixed_Subdivisions;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dynamic_mixed_subdivisions.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift