Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/dynamic_mixed_subdivisions.adb, Revision 1.1.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