Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Dynlift/triangulations.adb, Revision 1.1.1.1
1.1 maekawa 1: with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
2: with Standard_Random_Numbers; use Standard_Random_Numbers;
3:
4: --with text_io,Simplices_io; use text_io,Simplices_io;
5: --with Integer_Vectors_io; use Integer_Vectors_io;
6: --with integer_io; use integer_io;
7:
8: package body Triangulations is
9:
10: -- AUXILIAIRY HASH TABLE FOR Update_One :
11:
12: type Matrix_of_Triangulations
13: is array ( integer range <>, integer range <> ) of Triangulation;
14:
15: type Hash_Table ( n : natural ) is
16: record
17: weight1,weight2 : Link_to_Vector;
18: data : Matrix_of_Triangulations(0..n,0..n);
19: end record;
20:
21: procedure Init ( ht : in out Hash_Table; first,last : in integer ) is
22:
23: -- Initializes the hash table with Null_Triangulation.
24: -- Determines the weights for the hash function.
25: -- The numbers first and last determine the length of the weight function.
26:
27: begin
28: ht.weight1 := new Vector(first..last);
29: ht.weight2 := new Vector(first..last);
30: for i in first..last loop
31: ht.weight1(i) := Random(-last,last);
32: ht.weight2(i) := Random(-last,last);
33: end loop;
34: ht.weight1(last) := 1; -- TO AVOID PROBLEMS WITH HIGH LIFTING
35: ht.weight2(last) := 1;
36: for i in ht.data'range(1) loop
37: for j in ht.data'range(2) loop
38: ht.data(i,j) := Null_Triangulation;
39: end loop;
40: end loop;
41: end Init;
42:
43: procedure Hash_Function ( ht : in Hash_Table; s : in Simplex;
44: i,j : out integer ) is
45:
46: -- DESCRIPTION :
47: -- Computes the entries in the hash table for the simplex s.
48:
49: r1 : constant integer := ht.weight1.all*Vertex(s,1);
50: r2 : constant integer := ht.weight2.all*Vertex(s,2);
51:
52: begin
53: i := r1 mod (ht.n + 1);
54: j := r2 mod (ht.n + 1);
55: end Hash_Function;
56:
57: procedure Update ( ht : in out Hash_Table; s : in Simplex ) is
58:
59: -- DESCRIPTION :
60: -- Updates the hash table with the given simplex s.
61:
62: i,j : integer;
63:
64: begin
65: Hash_Function(ht,s,i,j);
66: Construct(s,ht.data(i,j));
67: end Update;
68:
69: function Is_In ( ht : Hash_Table; s : Simplex ) return boolean is
70:
71: -- DESCRIPTION :
72: -- Checks whether the simplex belongs to the hash table.
73:
74: i,j : integer;
75:
76: begin
77: Hash_Function(ht,s,i,j);
78: return Is_In(ht.data(i,j),s);
79: end Is_In;
80:
81: -- procedure Write_Distribution ( ht : in Hash_Table ) is
82: --
83: -- -- DESCRIPTION :
84: -- -- Writes the lengths of the lists in ht.
85: --
86: -- begin
87: -- for i in ht.data'range(1) loop
88: -- for j in ht.data'range(2) loop
89: -- put(Length_Of(ht.data(i,j)),1); put(' ');
90: -- end loop;
91: -- new_line;
92: -- end loop;
93: -- end Write_Distribution;
94:
95: procedure Clear ( ht : in out Hash_Table ) is
96:
97: -- DESCRIPTION :
98: -- Clears the allocated memory space for the hash table,
99: -- only a shallow copy is performed.
100:
101: begin
102: Clear(ht.weight1);
103: Clear(ht.weight2);
104: for i in ht.data'range(1) loop
105: for j in ht.data'range(2) loop
106: Lists_of_Simplices.Clear(Lists_of_Simplices.List(ht.data(i,j)));
107: end loop;
108: end loop;
109: end Clear;
110:
111: -- CREATORS :
112:
113: function Create ( s : Simplex ) return Triangulation is
114:
115: res : Triangulation;
116:
117: begin
118: Construct(s,res);
119: return res;
120: end Create;
121:
122: function Is_Inner ( n : natural; s : Simplex ) return boolean is
123: begin
124: for k in 1..n loop
125: if Neighbor(s,k) = Null_Simplex
126: then return false;
127: end if;
128: end loop;
129: return true;
130: end Is_Inner;
131:
132: procedure Update ( t : in out Triangulation; s : in out Simplex;
133: x : in Link_to_Vector ) is
134:
135: pos : vector(0..x'last) := Position(s,x.all);
136:
137: begin
138: for k in x'range loop
139: if Neighbor(s,k) = Null_Simplex
140: then if pos(k-1)*pos(pos'last) > 0
141: then Update(s,x,k);
142: Construct(Neighbor(s,k),t);
143: end if;
144: end if;
145: end loop;
146: end Update;
147:
148: procedure Connect_and_Update
149: ( s : in out Simplex; t : in out Triangulation ) is
150:
151: -- DESCRIPTION :
152: -- Connects the simplex with each other simplex in the triangulation
153: -- and adds it to the list t.
154:
155: tmp : Triangulation := t;
156:
157: begin
158: while not Is_Null(tmp) loop
159: declare
160: s2 : Simplex := Head_Of(tmp);
161: begin
162: Connect(s,s2);
163: --Set_Head(tmp,s2);
164: end;
165: tmp := Tail_Of(tmp);
166: end loop;
167: Construct(s,t);
168: end Connect_and_Update;
169:
170: procedure Connect_and_Update
171: ( t1 : in Triangulation; t2 : in out Triangulation ) is
172:
173: -- DESCRIPTION :
174: -- Connects all simplices in t1 properly with each other
175: -- and adds them to the list t2.
176:
177: tmp1 : Triangulation := t1;
178: tmp2 : Triangulation;
179: s1,s2 : Simplex;
180:
181: begin
182: while not Is_Null(tmp1) loop
183: s1 := Head_Of(tmp1);
184: tmp2 := Tail_Of(tmp1);
185: while not Is_Null(tmp2) loop
186: s2 := Head_Of(tmp2);
187: Connect(s1,s2);
188: tmp2 := Tail_Of(tmp2);
189: end loop;
190: tmp1 := Tail_Of(tmp1);
191: Construct(s1,t2);
192: end loop;
193: end Connect_and_Update;
194:
195: procedure Update ( t : in out Triangulation; x : in Link_to_Vector;
196: newt : out Triangulation ) is
197:
198: tmp : Triangulation := t;
199: s : Simplex;
200: nwt : Triangulation; -- for the new simplices that will contain x
201:
202: begin
203: while not Is_Null(tmp) loop
204: s := Head_Of(tmp);
205: if not Is_Inner(x'last,s)
206: then Update(nwt,s,x);
207: end if;
208: tmp := Tail_Of(tmp);
209: end loop;
210: Connect_and_Update(nwt,t);
211: newt := nwt;
212: end Update;
213:
214: procedure Update ( t : in out Triangulation; x : in Link_to_Vector ) is
215:
216: nwt : Triangulation; -- for the new simplices that will contain x
217:
218: begin
219: Update(t,x,nwt);
220: -- SHALLOW CLEAR :
221: Lists_of_Simplices.Clear(Lists_of_Simplices.List(nwt));
222: end Update;
223:
224: procedure Collect_Vertices
225: ( s : in Simplex; l : in out List ) is
226:
227: -- DESCRIPTION :
228: -- Collects the vertices in s and puts them in the list l.
229: -- Auxiliary routine for Update_One.
230:
231: pts : constant VecVec := Vertices(s);
232:
233: begin
234: for k in pts'range loop
235: if not Is_In(l,pts(k).all)
236: then declare
237: newpt : Link_to_Vector := new Vector'(pts(k).all);
238: begin
239: Construct(newpt,l);
240: end;
241: end if;
242: end loop;
243: end Collect_Vertices;
244:
245: function Has_Vertices_in_List ( s : Simplex; l : List ) return boolean is
246:
247: -- DESCRIPTION :
248: -- Returns true if the simplex s has vertices in the list l.
249:
250: pts : constant VecVec := Vertices(s);
251:
252: begin
253: for k in pts'range loop
254: if Is_In(l,pts(k).all)
255: then return true;
256: end if;
257: end loop;
258: return false;
259: end Has_Vertices_in_List;
260:
261: procedure Update_One
262: ( t : in out Triangulation; x : in Link_to_Vector ) is
263:
264: s : Simplex := Head_Of(t);
265: news,nei : Simplex;
266: nwt,leaves : Triangulation;
267: root : Hash_Table(2*x'last-1);
268: pos : Vector(0..x'last);
269: border_vertices : List;
270:
271: procedure Process_Tree_of_Updates is
272:
273: -- DESCRIPTION :
274: -- Constructs a tree of adjencencies which consists of
275: -- connected simplices on the edge of the triangulation.
276:
277: tmp,newleaves : Triangulation;
278: si,neisi : Simplex;
279:
280: begin
281: -- INITIALIZATION :
282: --Construct(s,root);
283: Init(root,x'first,x'last);
284: Update(root,s);
285: for k in x'range loop
286: declare
287: nei : Simplex := Neighbor(s,k);
288: begin
289: if (nei /= Null_Simplex) and then not Is_Vertex(nei,x.all)
290: then Construct(nei,leaves);
291: end if;
292: end;
293: end loop;
294: -- PERFORM THE UPDATES :
295: loop
296: -- INVARIANT CONDITION :
297: -- the simplices considered have vertices on the border close to x
298: newleaves := Null_Triangulation;
299: tmp := leaves;
300: while not Is_Null(tmp) loop
301: si := Head_Of(tmp);
302: -- UPDATE root TO PREVENT TURNING BACK :
303: Update(root,si);
304: -- COMPUTE NEW SIMPLICES AND LEAVES :
305: if not Is_Inner(x'last,si)
306: then pos := Position(si,x.all);
307: for k in x'range loop
308: -- LOOK IN THE DIRECTION TOWARDS x :
309: if pos(k-1)*pos(pos'last) > 0
310: then neisi := Neighbor(si,k);
311: if neisi = Null_Simplex
312: then Update(si,x,k); -- NEW SIMPLEX
313: neisi := Neighbor(si,k);
314: Construct(neisi,nwt);
315: Collect_Vertices(neisi,border_vertices);
316: end if;
317: end if;
318: end loop;
319: end if;
320: for l in x'range loop -- GO FURTHER
321: neisi := Neighbor(si,l);
322: if (neisi /= Null_Simplex)
323: and then not Is_Vertex(neisi,x.all)
324: and then Has_Vertices_in_List(neisi,border_vertices)
325: and then not Is_In(leaves,neisi)
326: and then not Is_In(newleaves,neisi)
327: and then not Is_In(root,neisi)
328: then Construct(neisi,newleaves);
329: end if;
330: end loop;
331: tmp := Tail_Of(tmp);
332: end loop;
333: exit when (newleaves = Null_Triangulation);
334: Lists_of_Simplices.Clear(Lists_of_Simplices.List(leaves));
335: leaves := newleaves;
336: end loop;
337: -- SHALLOW CLEAR :
338: -- Lists_of_Simplices.Clear(Lists_of_Simplices.List(root));
339: -- put_line("Distribution of root : "); Write_Distribution(root);
340: Clear(root);
341: Lists_of_Simplices.Clear(Lists_of_Simplices.List(leaves));
342: end Process_Tree_of_Updates;
343:
344: begin
345: nei := s; -- USED TO SEE WHETHER s WILL CHANGE
346: pos := Position(s,x.all);
347: Update_One(s,x,pos,news);
348: Construct(news,nwt);
349: Collect_Vertices(news,border_vertices);
350: if s /= nei
351: then pos := Position(s,x.all); -- BECAUSE s IS CHANGED !!!
352: end if;
353: for k in x'range loop -- COMPUTE OTHER NEIGHBORS
354: if pos(k-1)*pos(pos'last) > 0
355: then nei := Neighbor(s,k);
356: if nei = Null_Simplex
357: then Update(s,x,k);
358: nei := Neighbor(s,k);
359: Construct(nei,nwt);
360: Collect_Vertices(nei,border_vertices);
361: end if;
362: end if;
363: end loop;
364: Process_Tree_of_Updates;
365: Clear(border_vertices);
366: Connect_and_Update(nwt,t);
367: Lists_of_Simplices.Clear(Lists_of_Simplices.List(nwt));
368: end Update_One;
369:
370: procedure Connect ( t : in out Triangulation ) is
371:
372: tmp1 : Triangulation := t;
373: tmp2 : Triangulation;
374: s1,s2 : Simplex;
375:
376: begin
377: while not Is_Null(tmp1) loop
378: s1 := Head_Of(tmp1);
379: tmp2 := Tail_Of(tmp1);
380: while not Is_Null(tmp2) loop
381: s2 := Head_Of(tmp2);
382: Connect(s1,s2);
383: --Set_Head(tmp2,s2);
384: tmp2 := Tail_Of(tmp2);
385: end loop;
386: --Set_Head(tmp1,s1);
387: tmp1 := Tail_Of(tmp1);
388: end loop;
389: end Connect;
390:
391: -- THE OPTIMAL DISCRETE CONSERVATIVE LIFTING FUNCTION :
392:
393: -- procedure Write_Neighbors ( s : in Simplex ) is
394: -- begin
395: -- put("The normal of simplex s : "); put(Normal(s)); new_line;
396: -- put_line(" with vertices : "); put(s);
397: -- for k in Normal(s)'range loop
398: -- if Neighbor(s,k) /= Null_Simplex
399: -- then put("normal of a not null neighbors of s :");
400: -- put(Normal(Neighbor(s,k))); new_line;
401: -- end if;
402: -- end loop;
403: -- end Write_Neighbors;
404:
405: -- procedure Write ( t : in Triangulation ) is
406: --
407: -- tmp : Triangulation := t;
408: --
409: -- begin
410: -- while not Is_Null(tmp) loop
411: -- Write_Neighbors(Head_Of(tmp));
412: -- tmp := Tail_Of(tmp);
413: -- end loop;
414: -- end Write;
415:
416: procedure Flatten ( t : in out Triangulation ) is
417:
418: tmp : Triangulation := t;
419: s : Simplex;
420:
421: -- IMPLEMENTATION REQUIREMENT :
422: -- Cells that are already flattened are grouped at the end of the list.
423:
424: begin
425: -- put_line("THE TRIANGULATION BEFORE LIFTING : "); Write(t);
426: while not Is_Null(tmp) loop
427: s := Head_Of(tmp);
428: exit when Is_Flat(s);
429: -- put_line("NEIGHBORS BEFORE FLATTENING : "); Write_Neighbors(s);
430: Flatten(s);
431: -- put_line("NEIGHBORS AFTER FLATTENING : "); Write_Neighbors(s);
432: --Set_Head(tmp,s);
433: tmp := Tail_Of(tmp);
434: end loop;
435: -- put_line("THE TRIANGULATION AFTER LIFTING : "); Write(t);
436: end Flatten;
437:
438: -- SELECTORS :
439:
440: function Is_Vertex ( t : Triangulation; x : Vector ) return boolean is
441:
442: tmp : Triangulation := t;
443:
444: begin
445: while not Is_Null(tmp) loop
446: if Is_Vertex(Head_Of(tmp),x)
447: then return true;
448: end if;
449: tmp := Tail_Of(tmp);
450: end loop;
451: return false;
452: end Is_Vertex;
453:
454: function Vertices ( t : Triangulation ) return List is
455:
456: -- DESCRIPTION :
457: -- Returns a list with all vertices in the simplices of t.
458:
459: res,res_last : List;
460: tmp : Triangulation := t;
461:
462: begin
463: res_last := res;
464: while not Is_Null(tmp) loop
465: declare
466: s : Simplex := Head_Of(tmp);
467: v : constant VecVec := Vertices(s);
468: begin
469: for i in v'range loop
470: if not Is_In(res,v(i))
471: then Append(res,res_last,v(i));
472: end if;
473: end loop;
474: end;
475: tmp := Tail_Of(tmp);
476: end loop;
477: return res;
478: end Vertices;
479:
480: function Vertices ( t : Triangulation; x : vector ) return List is
481:
482: res,res_last : List;
483: tmp : Triangulation := t;
484:
485: begin
486: res_last := res;
487: while not Is_Null(tmp) loop
488: declare
489: s : Simplex := Head_Of(tmp);
490: begin
491: if Is_Vertex(s,x)
492: then
493: declare
494: v : constant VecVec := Vertices(s);
495: begin
496: for i in v'range loop
497: if not Is_In(res,v(i))
498: then Append(res,res_last,v(i));
499: end if;
500: end loop;
501: end;
502: end if;
503: tmp := Tail_Of(tmp);
504: end;
505: end loop;
506: return res;
507: end Vertices;
508:
509: function Vertices ( t : Triangulation ) return VecVec is
510:
511: vertri : List := Vertices(t);
512: res : VecVec(1..Length_Of(vertri));
513: tmp : List := vertri;
514:
515: begin
516: for i in res'range loop
517: res(i) := Head_Of(tmp);
518: tmp := Tail_Of(tmp);
519: end loop;
520: Shallow_Clear(vertri);
521: -- Lists_of_Link_to_Integer_Vectors.Clear
522: -- (Lists_of_Link_to_Integer_Vectors.List(vertri));
523: return res;
524: end Vertices;
525:
526: function Vertices ( t : Triangulation; x : vector ) return VecVec is
527:
528: vertri : List := Vertices(t,x);
529: res : VecVec(1..Length_Of(vertri));
530: tmp : List := vertri;
531:
532: begin
533: for i in res'range loop
534: res(i) := Head_Of(tmp);
535: tmp := Tail_Of(tmp);
536: end loop;
537: Shallow_Clear(vertri);
538: -- Lists_of_Link_to_Integer_Vectors.Clear
539: -- (Lists_of_Link_to_Integer_Vectors.List(vertri));
540: return res;
541: end Vertices;
542:
543: function Is_In1 ( t : Triangulation; x : Vector ) return boolean is
544:
545: tmp : Triangulation := t;
546:
547: begin
548: while not Is_Null(tmp) loop
549: if Is_In(Head_Of(tmp),x)
550: then return true;
551: end if;
552: tmp := Tail_Of(tmp);
553: end loop;
554: return false;
555: end Is_In1;
556:
557: function Is_In ( t : Triangulation; x : Vector ) return boolean is
558: begin
559: return Is_In_All(Head_Of(t),x);
560: end Is_In;
561:
562: function Is_In ( t : Triangulation; x : Vector ) return Simplex is
563: begin
564: return Is_In_All(Head_Of(t),x);
565: end Is_In;
566:
567: function Is_In ( t : Triangulation; s : Simplex ) return boolean is
568:
569: tmp : Triangulation := t;
570:
571: begin
572: -- put("Length of triangulation : "); put(Length_Of(t),1); new_line;
573: while not Is_Null(tmp) loop
574: if Equal(s,Head_Of(tmp))
575: then return true;
576: else tmp := Tail_Of(tmp);
577: end if;
578: end loop;
579: return false;
580: end Is_In;
581:
582: function Volume ( t : Triangulation ) return natural is
583:
584: tmp : Triangulation := t;
585: vol : natural := 0;
586:
587: begin
588: while not Is_Null(tmp) loop
589: vol := vol + Volume(Head_Of(tmp));
590: tmp := Tail_Of(tmp);
591: end loop;
592: return vol;
593: end Volume;
594:
595: -- DESTRUCTOR :
596:
597: procedure Clear ( t : in out Triangulation ) is
598:
599: tmp : Triangulation := t;
600:
601: begin
602: while not Is_Null(tmp) loop
603: declare
604: s : Simplex := Head_Of(tmp);
605: begin
606: Clear(s);
607: end;
608: tmp := Tail_Of(tmp);
609: end loop;
610: Lists_of_Simplices.Clear(Lists_of_Simplices.List(t));
611: end Clear;
612:
613: end Triangulations;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to triangulations.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