File: [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift / triangulations.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:28 2000 UTC (23 years, 10 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with Lists_of_Integer_Vectors; use Lists_of_Integer_Vectors;
with Standard_Random_Numbers; use Standard_Random_Numbers;
--with text_io,Simplices_io; use text_io,Simplices_io;
--with Integer_Vectors_io; use Integer_Vectors_io;
--with integer_io; use integer_io;
package body Triangulations is
-- AUXILIAIRY HASH TABLE FOR Update_One :
type Matrix_of_Triangulations
is array ( integer range <>, integer range <> ) of Triangulation;
type Hash_Table ( n : natural ) is
record
weight1,weight2 : Link_to_Vector;
data : Matrix_of_Triangulations(0..n,0..n);
end record;
procedure Init ( ht : in out Hash_Table; first,last : in integer ) is
-- Initializes the hash table with Null_Triangulation.
-- Determines the weights for the hash function.
-- The numbers first and last determine the length of the weight function.
begin
ht.weight1 := new Vector(first..last);
ht.weight2 := new Vector(first..last);
for i in first..last loop
ht.weight1(i) := Random(-last,last);
ht.weight2(i) := Random(-last,last);
end loop;
ht.weight1(last) := 1; -- TO AVOID PROBLEMS WITH HIGH LIFTING
ht.weight2(last) := 1;
for i in ht.data'range(1) loop
for j in ht.data'range(2) loop
ht.data(i,j) := Null_Triangulation;
end loop;
end loop;
end Init;
procedure Hash_Function ( ht : in Hash_Table; s : in Simplex;
i,j : out integer ) is
-- DESCRIPTION :
-- Computes the entries in the hash table for the simplex s.
r1 : constant integer := ht.weight1.all*Vertex(s,1);
r2 : constant integer := ht.weight2.all*Vertex(s,2);
begin
i := r1 mod (ht.n + 1);
j := r2 mod (ht.n + 1);
end Hash_Function;
procedure Update ( ht : in out Hash_Table; s : in Simplex ) is
-- DESCRIPTION :
-- Updates the hash table with the given simplex s.
i,j : integer;
begin
Hash_Function(ht,s,i,j);
Construct(s,ht.data(i,j));
end Update;
function Is_In ( ht : Hash_Table; s : Simplex ) return boolean is
-- DESCRIPTION :
-- Checks whether the simplex belongs to the hash table.
i,j : integer;
begin
Hash_Function(ht,s,i,j);
return Is_In(ht.data(i,j),s);
end Is_In;
-- procedure Write_Distribution ( ht : in Hash_Table ) is
--
-- -- DESCRIPTION :
-- -- Writes the lengths of the lists in ht.
--
-- begin
-- for i in ht.data'range(1) loop
-- for j in ht.data'range(2) loop
-- put(Length_Of(ht.data(i,j)),1); put(' ');
-- end loop;
-- new_line;
-- end loop;
-- end Write_Distribution;
procedure Clear ( ht : in out Hash_Table ) is
-- DESCRIPTION :
-- Clears the allocated memory space for the hash table,
-- only a shallow copy is performed.
begin
Clear(ht.weight1);
Clear(ht.weight2);
for i in ht.data'range(1) loop
for j in ht.data'range(2) loop
Lists_of_Simplices.Clear(Lists_of_Simplices.List(ht.data(i,j)));
end loop;
end loop;
end Clear;
-- CREATORS :
function Create ( s : Simplex ) return Triangulation is
res : Triangulation;
begin
Construct(s,res);
return res;
end Create;
function Is_Inner ( n : natural; s : Simplex ) return boolean is
begin
for k in 1..n loop
if Neighbor(s,k) = Null_Simplex
then return false;
end if;
end loop;
return true;
end Is_Inner;
procedure Update ( t : in out Triangulation; s : in out Simplex;
x : in Link_to_Vector ) is
pos : vector(0..x'last) := Position(s,x.all);
begin
for k in x'range loop
if Neighbor(s,k) = Null_Simplex
then if pos(k-1)*pos(pos'last) > 0
then Update(s,x,k);
Construct(Neighbor(s,k),t);
end if;
end if;
end loop;
end Update;
procedure Connect_and_Update
( s : in out Simplex; t : in out Triangulation ) is
-- DESCRIPTION :
-- Connects the simplex with each other simplex in the triangulation
-- and adds it to the list t.
tmp : Triangulation := t;
begin
while not Is_Null(tmp) loop
declare
s2 : Simplex := Head_Of(tmp);
begin
Connect(s,s2);
--Set_Head(tmp,s2);
end;
tmp := Tail_Of(tmp);
end loop;
Construct(s,t);
end Connect_and_Update;
procedure Connect_and_Update
( t1 : in Triangulation; t2 : in out Triangulation ) is
-- DESCRIPTION :
-- Connects all simplices in t1 properly with each other
-- and adds them to the list t2.
tmp1 : Triangulation := t1;
tmp2 : Triangulation;
s1,s2 : Simplex;
begin
while not Is_Null(tmp1) loop
s1 := Head_Of(tmp1);
tmp2 := Tail_Of(tmp1);
while not Is_Null(tmp2) loop
s2 := Head_Of(tmp2);
Connect(s1,s2);
tmp2 := Tail_Of(tmp2);
end loop;
tmp1 := Tail_Of(tmp1);
Construct(s1,t2);
end loop;
end Connect_and_Update;
procedure Update ( t : in out Triangulation; x : in Link_to_Vector;
newt : out Triangulation ) is
tmp : Triangulation := t;
s : Simplex;
nwt : Triangulation; -- for the new simplices that will contain x
begin
while not Is_Null(tmp) loop
s := Head_Of(tmp);
if not Is_Inner(x'last,s)
then Update(nwt,s,x);
end if;
tmp := Tail_Of(tmp);
end loop;
Connect_and_Update(nwt,t);
newt := nwt;
end Update;
procedure Update ( t : in out Triangulation; x : in Link_to_Vector ) is
nwt : Triangulation; -- for the new simplices that will contain x
begin
Update(t,x,nwt);
-- SHALLOW CLEAR :
Lists_of_Simplices.Clear(Lists_of_Simplices.List(nwt));
end Update;
procedure Collect_Vertices
( s : in Simplex; l : in out List ) is
-- DESCRIPTION :
-- Collects the vertices in s and puts them in the list l.
-- Auxiliary routine for Update_One.
pts : constant VecVec := Vertices(s);
begin
for k in pts'range loop
if not Is_In(l,pts(k).all)
then declare
newpt : Link_to_Vector := new Vector'(pts(k).all);
begin
Construct(newpt,l);
end;
end if;
end loop;
end Collect_Vertices;
function Has_Vertices_in_List ( s : Simplex; l : List ) return boolean is
-- DESCRIPTION :
-- Returns true if the simplex s has vertices in the list l.
pts : constant VecVec := Vertices(s);
begin
for k in pts'range loop
if Is_In(l,pts(k).all)
then return true;
end if;
end loop;
return false;
end Has_Vertices_in_List;
procedure Update_One
( t : in out Triangulation; x : in Link_to_Vector ) is
s : Simplex := Head_Of(t);
news,nei : Simplex;
nwt,leaves : Triangulation;
root : Hash_Table(2*x'last-1);
pos : Vector(0..x'last);
border_vertices : List;
procedure Process_Tree_of_Updates is
-- DESCRIPTION :
-- Constructs a tree of adjencencies which consists of
-- connected simplices on the edge of the triangulation.
tmp,newleaves : Triangulation;
si,neisi : Simplex;
begin
-- INITIALIZATION :
--Construct(s,root);
Init(root,x'first,x'last);
Update(root,s);
for k in x'range loop
declare
nei : Simplex := Neighbor(s,k);
begin
if (nei /= Null_Simplex) and then not Is_Vertex(nei,x.all)
then Construct(nei,leaves);
end if;
end;
end loop;
-- PERFORM THE UPDATES :
loop
-- INVARIANT CONDITION :
-- the simplices considered have vertices on the border close to x
newleaves := Null_Triangulation;
tmp := leaves;
while not Is_Null(tmp) loop
si := Head_Of(tmp);
-- UPDATE root TO PREVENT TURNING BACK :
Update(root,si);
-- COMPUTE NEW SIMPLICES AND LEAVES :
if not Is_Inner(x'last,si)
then pos := Position(si,x.all);
for k in x'range loop
-- LOOK IN THE DIRECTION TOWARDS x :
if pos(k-1)*pos(pos'last) > 0
then neisi := Neighbor(si,k);
if neisi = Null_Simplex
then Update(si,x,k); -- NEW SIMPLEX
neisi := Neighbor(si,k);
Construct(neisi,nwt);
Collect_Vertices(neisi,border_vertices);
end if;
end if;
end loop;
end if;
for l in x'range loop -- GO FURTHER
neisi := Neighbor(si,l);
if (neisi /= Null_Simplex)
and then not Is_Vertex(neisi,x.all)
and then Has_Vertices_in_List(neisi,border_vertices)
and then not Is_In(leaves,neisi)
and then not Is_In(newleaves,neisi)
and then not Is_In(root,neisi)
then Construct(neisi,newleaves);
end if;
end loop;
tmp := Tail_Of(tmp);
end loop;
exit when (newleaves = Null_Triangulation);
Lists_of_Simplices.Clear(Lists_of_Simplices.List(leaves));
leaves := newleaves;
end loop;
-- SHALLOW CLEAR :
-- Lists_of_Simplices.Clear(Lists_of_Simplices.List(root));
-- put_line("Distribution of root : "); Write_Distribution(root);
Clear(root);
Lists_of_Simplices.Clear(Lists_of_Simplices.List(leaves));
end Process_Tree_of_Updates;
begin
nei := s; -- USED TO SEE WHETHER s WILL CHANGE
pos := Position(s,x.all);
Update_One(s,x,pos,news);
Construct(news,nwt);
Collect_Vertices(news,border_vertices);
if s /= nei
then pos := Position(s,x.all); -- BECAUSE s IS CHANGED !!!
end if;
for k in x'range loop -- COMPUTE OTHER NEIGHBORS
if pos(k-1)*pos(pos'last) > 0
then nei := Neighbor(s,k);
if nei = Null_Simplex
then Update(s,x,k);
nei := Neighbor(s,k);
Construct(nei,nwt);
Collect_Vertices(nei,border_vertices);
end if;
end if;
end loop;
Process_Tree_of_Updates;
Clear(border_vertices);
Connect_and_Update(nwt,t);
Lists_of_Simplices.Clear(Lists_of_Simplices.List(nwt));
end Update_One;
procedure Connect ( t : in out Triangulation ) is
tmp1 : Triangulation := t;
tmp2 : Triangulation;
s1,s2 : Simplex;
begin
while not Is_Null(tmp1) loop
s1 := Head_Of(tmp1);
tmp2 := Tail_Of(tmp1);
while not Is_Null(tmp2) loop
s2 := Head_Of(tmp2);
Connect(s1,s2);
--Set_Head(tmp2,s2);
tmp2 := Tail_Of(tmp2);
end loop;
--Set_Head(tmp1,s1);
tmp1 := Tail_Of(tmp1);
end loop;
end Connect;
-- THE OPTIMAL DISCRETE CONSERVATIVE LIFTING FUNCTION :
-- procedure Write_Neighbors ( s : in Simplex ) is
-- begin
-- put("The normal of simplex s : "); put(Normal(s)); new_line;
-- put_line(" with vertices : "); put(s);
-- for k in Normal(s)'range loop
-- if Neighbor(s,k) /= Null_Simplex
-- then put("normal of a not null neighbors of s :");
-- put(Normal(Neighbor(s,k))); new_line;
-- end if;
-- end loop;
-- end Write_Neighbors;
-- procedure Write ( t : in Triangulation ) is
--
-- tmp : Triangulation := t;
--
-- begin
-- while not Is_Null(tmp) loop
-- Write_Neighbors(Head_Of(tmp));
-- tmp := Tail_Of(tmp);
-- end loop;
-- end Write;
procedure Flatten ( t : in out Triangulation ) is
tmp : Triangulation := t;
s : Simplex;
-- IMPLEMENTATION REQUIREMENT :
-- Cells that are already flattened are grouped at the end of the list.
begin
-- put_line("THE TRIANGULATION BEFORE LIFTING : "); Write(t);
while not Is_Null(tmp) loop
s := Head_Of(tmp);
exit when Is_Flat(s);
-- put_line("NEIGHBORS BEFORE FLATTENING : "); Write_Neighbors(s);
Flatten(s);
-- put_line("NEIGHBORS AFTER FLATTENING : "); Write_Neighbors(s);
--Set_Head(tmp,s);
tmp := Tail_Of(tmp);
end loop;
-- put_line("THE TRIANGULATION AFTER LIFTING : "); Write(t);
end Flatten;
-- SELECTORS :
function Is_Vertex ( t : Triangulation; x : Vector ) return boolean is
tmp : Triangulation := t;
begin
while not Is_Null(tmp) loop
if Is_Vertex(Head_Of(tmp),x)
then return true;
end if;
tmp := Tail_Of(tmp);
end loop;
return false;
end Is_Vertex;
function Vertices ( t : Triangulation ) return List is
-- DESCRIPTION :
-- Returns a list with all vertices in the simplices of t.
res,res_last : List;
tmp : Triangulation := t;
begin
res_last := res;
while not Is_Null(tmp) loop
declare
s : Simplex := Head_Of(tmp);
v : constant VecVec := Vertices(s);
begin
for i in v'range loop
if not Is_In(res,v(i))
then Append(res,res_last,v(i));
end if;
end loop;
end;
tmp := Tail_Of(tmp);
end loop;
return res;
end Vertices;
function Vertices ( t : Triangulation; x : vector ) return List is
res,res_last : List;
tmp : Triangulation := t;
begin
res_last := res;
while not Is_Null(tmp) loop
declare
s : Simplex := Head_Of(tmp);
begin
if Is_Vertex(s,x)
then
declare
v : constant VecVec := Vertices(s);
begin
for i in v'range loop
if not Is_In(res,v(i))
then Append(res,res_last,v(i));
end if;
end loop;
end;
end if;
tmp := Tail_Of(tmp);
end;
end loop;
return res;
end Vertices;
function Vertices ( t : Triangulation ) return VecVec is
vertri : List := Vertices(t);
res : VecVec(1..Length_Of(vertri));
tmp : List := vertri;
begin
for i in res'range loop
res(i) := Head_Of(tmp);
tmp := Tail_Of(tmp);
end loop;
Shallow_Clear(vertri);
-- Lists_of_Link_to_Integer_Vectors.Clear
-- (Lists_of_Link_to_Integer_Vectors.List(vertri));
return res;
end Vertices;
function Vertices ( t : Triangulation; x : vector ) return VecVec is
vertri : List := Vertices(t,x);
res : VecVec(1..Length_Of(vertri));
tmp : List := vertri;
begin
for i in res'range loop
res(i) := Head_Of(tmp);
tmp := Tail_Of(tmp);
end loop;
Shallow_Clear(vertri);
-- Lists_of_Link_to_Integer_Vectors.Clear
-- (Lists_of_Link_to_Integer_Vectors.List(vertri));
return res;
end Vertices;
function Is_In1 ( t : Triangulation; x : Vector ) return boolean is
tmp : Triangulation := t;
begin
while not Is_Null(tmp) loop
if Is_In(Head_Of(tmp),x)
then return true;
end if;
tmp := Tail_Of(tmp);
end loop;
return false;
end Is_In1;
function Is_In ( t : Triangulation; x : Vector ) return boolean is
begin
return Is_In_All(Head_Of(t),x);
end Is_In;
function Is_In ( t : Triangulation; x : Vector ) return Simplex is
begin
return Is_In_All(Head_Of(t),x);
end Is_In;
function Is_In ( t : Triangulation; s : Simplex ) return boolean is
tmp : Triangulation := t;
begin
-- put("Length of triangulation : "); put(Length_Of(t),1); new_line;
while not Is_Null(tmp) loop
if Equal(s,Head_Of(tmp))
then return true;
else tmp := Tail_Of(tmp);
end if;
end loop;
return false;
end Is_In;
function Volume ( t : Triangulation ) return natural is
tmp : Triangulation := t;
vol : natural := 0;
begin
while not Is_Null(tmp) loop
vol := vol + Volume(Head_Of(tmp));
tmp := Tail_Of(tmp);
end loop;
return vol;
end Volume;
-- DESTRUCTOR :
procedure Clear ( t : in out Triangulation ) is
tmp : Triangulation := t;
begin
while not Is_Null(tmp) loop
declare
s : Simplex := Head_Of(tmp);
begin
Clear(s);
end;
tmp := Tail_Of(tmp);
end loop;
Lists_of_Simplices.Clear(Lists_of_Simplices.List(t));
end Clear;
end Triangulations;