File: [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Dynlift / dynamic_mixed_subdivisions.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 Standard_Integer_VecVecs; use Standard_Integer_VecVecs;
with Transforming_Integer_Vector_Lists; use Transforming_Integer_Vector_Lists;
with Initial_Mixed_Cell;
with Vertices,Simplices; use Vertices,Simplices;
with Global_Dynamic_Triangulation; use Global_Dynamic_Triangulation;
with Integer_Lifting_Utilities; use Integer_Lifting_Utilities;
with Dynamic_Lifting_Functions; use Dynamic_Lifting_Functions;
with Flatten_Mixed_Subdivisions; use Flatten_Mixed_Subdivisions;
with Enumerate_Faces_of_Polytope; use Enumerate_Faces_of_Polytope;
with Common_Faces_of_Polytope; use Common_Faces_of_Polytope;
with Integer_Pruning_Methods; use Integer_Pruning_Methods;
with Mixed_Volume_Computation; use Mixed_Volume_Computation;
with Contributions_to_Mixed_Volume; use Contributions_to_Mixed_Volume;
package body Dynamic_Mixed_Subdivisions is
-- UTILITIES :
function Is_Empty ( points : Array_of_Lists ) return boolean is
begin
for i in points'range loop
if not Is_Null(points(i))
then return false;
end if;
end loop;
return true;
end Is_Empty;
function First_Non_Empty ( points : Array_of_Lists ) return integer is
-- DESCRIPTION :
-- Returns the index of the first non empty set in the points.
res : integer := points'first - 1;
begin
for i in points'range loop
if not Is_Null(points(i))
then res := i;
end if;
exit when res >= points'first;
end loop;
return res;
end First_Non_Empty;
function Is_In ( fs : Face_Structure; point : vector ) return boolean is
begin
if Is_Null(fs.t)
then return Is_In(fs.l,point);
else return Is_In(fs.t,point);
end if;
end Is_In;
-- procedure put ( n : in natural; fs : in Face_Structure ) is
-- begin
-- put_line("The list of lifted points : "); put(fs.l);
-- put_line("The list of k-faces : "); put(fs.f);
-- put_line("The triangulation : "); put(n,fs.t);
-- end put;
-- procedure put ( n : in natural; fs : in Face_Structures ) is
-- begin
-- for i in fs'range loop
-- put("face structure for component "); put(i,1);
-- put_line(" :");
-- put(n,fs(i));
-- end loop;
-- end put;
-- FLATTENING :
procedure Flatten ( points : in out VecVec ) is
pt : Link_to_Vector;
begin
for i in points'range loop
pt := points(i);
if pt(pt'last) /= 0
then pt(pt'last) := 0;
end if;
end loop;
end Flatten;
procedure Flatten ( f : in out Face ) is
begin
Flatten(f.all);
end Flatten;
procedure Flatten ( fs : in out Faces ) is
tmp : Faces := fs;
f : Face;
begin
while not Is_Null(tmp) loop
f := Head_Of(tmp);
Flatten(f);
Set_Head(tmp,f);
tmp := Tail_Of(tmp);
end loop;
end Flatten;
procedure Flatten ( fs : in out Face_Structure ) is
begin
Flatten(fs.l); Flatten(fs.t); Flatten(fs.f);
end Flatten;
procedure Flatten ( fs : in out Face_Structures ) is
begin
for i in fs'range loop
Flatten(fs(i));
end loop;
end Flatten;
-- ZERO CONTRIBUTION :
function Extract_Facet ( n : natural; s : Simplex ) return Face is
ver : constant VecVec := Simplices.Vertices(s);
res : Face := new VecVec(ver'range);
begin
for i in res'range loop
res(i) := new Standard_Integer_Vectors.Vector'(ver(i)(1..n));
end loop;
return res;
end Extract_Facet;
function Extract_Facets ( n : natural; t : Triangulation; x : Vector )
return Faces is
-- DESCRIPTION :
-- Returns the list of facets in t that all contain x.
-- The facets are given in their original coordinates.
-- REQUIRED : x is lifted, i.e., x'length = n+1.
tmp : Triangulation := t;
res,res_last : Faces;
begin
while not Is_Null(tmp) loop
declare
s : Simplex := Head_Of(tmp);
begin
if Is_Vertex(s,x)
then Append(res,res_last,Extract_Facet(n,s));
end if;
end;
tmp := Tail_Of(tmp);
end loop;
return res;
end Extract_Facets;
function Zero_Contribution ( n : natural; fs : Face_Structures;
x : Vector; i : natural ) return boolean is
-- DESCRIPTION :
-- Returns true if the point x does not contribute to the mixed volume
-- when considered to the ith component of the already lifted points.
-- REQUIRED : x'length = n+1, n = dimension before lifting.
res : boolean := false;
supp : Array_of_Lists(fs'range);
begin
for j in supp'range loop
supp(j) := Reduce(fs(j).l,n+1);
end loop;
if Length_Of(fs(i).t) > 1
then declare
facets : Faces := Extract_Facets(n,fs(i).t,x);
begin
res := Simple_Zero_Contribution(supp,facets,x(1..n),i);
end;
else res := Simple_Zero_Contribution(supp,x(1..n),i);
end if;
Deep_Clear(supp);
return res;
end Zero_Contribution;
-- INITIALIZATION :
procedure Initialize
( n : in natural; mix : in Vector; points : in Array_of_Lists;
mixsub : in out Mixed_Subdivision;
fs : in out Face_Structures; rest : in out Array_of_Lists;
done : in out boolean ) is
-- DESCRIPTION :
-- Performs initialization of the dynamic lifting algorithm.
-- ON ENTRY :
-- n length of the vectors in points;
-- mix type of mixture;
-- mixsub mixed subdivision of the lifted points;
-- fs face structures of the lifted points.
-- ON RETURN :
-- mixsub if empty on entry, then it contains the initial mixed
-- cell, in case the problem is nondegenerate;
-- rest rest of point list to be processed;
-- done if true, then either the mixed volume equals zero,
-- or there are no more points to process.
null_points : Array_of_Lists(points'range);
begin
if Is_Null(mixsub)
then -- compute an initial mixed cell
declare
mic : Mixed_Cell;
begin
Initial_Mixed_Cell(n,mix,points,mic,rest);
if Mixed_Volume(n,mix,mic) > 0 -- check on degeneracy
then
-- initialize the mixed subdivision :
Construct(mic,mixsub);
-- initialize the face structures :
for i in mic.pts'range loop
declare
tmp : List := mic.pts(i);
fc : Face;
begin
while not Is_Null(tmp) loop
Append(fs(i).l,fs(i).last,Head_Of(tmp).all);
tmp := Tail_of(tmp);
end loop;
fc := new VecVec'(Shallow_Create(mic.pts(i)));
Construct(fc,fs(i).f);
end;
end loop;
if mix'last = mix'first
then declare
v : constant VecVec := Shallow_Create(fs(fs'first).l);
s : Simplex := Create(v);
begin
fs(fs'first).t := Create(s);
end;
end if;
done := Is_Empty(rest);
else -- degenerate problem => mixed volume = 0
rest := null_points; done := true;
end if;
end;
else
rest := points;
done := Is_Empty(rest);
end if;
end Initialize;
function Initial_Triangulation
( n : in natural; l : in List; point : in Link_to_Vector )
return Triangulation is
-- DESCRIPTION :
-- Given a list of lifted points with Length_Of(l) >= n and another
-- lifted point, a nonempty triangulation will be returned, when the
-- volume of \conv(l U {point}) > 0.
res : Triangulation;
del,span : List;
delpoint : Link_to_Vector := new vector(1..n);
function Deflate ( lifted : in List ) return List is
-- DESCRIPTION :
-- Discards the lifting value of the points in the list l.
res,res_last : List;
tmp : List := lifted;
pt : Link_to_Vector;
begin
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
declare
dpt : Link_to_Vector := new vector(1..n);
begin
dpt.all := pt(dpt'range);
Append(res,res_last,dpt);
end;
tmp := Tail_Of(tmp);
end loop;
return res;
end Deflate;
function Lift ( pt,liftpt : Link_to_Vector; lifted : List )
return Link_to_Vector is
-- DESCRIPTION :
-- Compares the value of pt with that of liftpt and
-- looks up the lifting value of the point pt in the list lifted.
-- The lifted point will be returned.
-- REQUIRED :
-- Either the point pt should be equal to Deflate(liftpt)
-- or it should belong to Deflate(lifted).
tmp : List := lifted;
res : Link_to_Vector;
begin
if pt.all = liftpt(pt'range)
then return liftpt;
else while not Is_Null(tmp) loop
res := Head_Of(tmp);
if pt.all = res(pt'range)
then return res;
else tmp := Tail_Of(tmp);
end if;
end loop;
return res;
end if;
end Lift;
begin
del := Deflate(l);
delpoint.all := point(delpoint'range);
Construct(delpoint,del);
span := Extremal_Points(n,del);
if Length_Of(span) = n+1
then
declare
verpts : VecVec(1..n+1);
tmp : List := span;
pt : Link_to_Vector;
s : Simplex;
begin
for i in verpts'range loop
verpts(i) := Lift(Head_Of(tmp),point,l);
tmp := Tail_Of(tmp);
end loop;
s := Create(verpts);
res := Create(s);
tmp := l;
while not Is_Null(tmp) loop
pt := Head_Of(tmp);
if not Is_In(span,pt)
then Update(res,pt);
end if;
tmp := Tail_Of(tmp);
end loop;
end;
end if;
Clear(del); Clear(span); Clear(delpoint);
return res;
end Initial_Triangulation;
-- CHOOSE NEXT POINT :
procedure Next_Point
( n : in natural; points : in out Array_of_Lists;
order : in boolean;
index : in out integer; point : out Link_to_Vector ) is
-- DESCRIPTION :
-- Chooses a point out of the lists of points.
-- ON ENTRY :
-- n length of the elements in the point lists;
-- points nonempty lists of points;
-- order if true, then the first element out of the next first
-- nonempty list will be chosen;
-- index indicates component where not Is_Null(points(index)).
-- ON RETURN :
-- points lists of points with the point removed;
-- index points to the next nonempty list,
-- if index < points'first, then Is_Empty(points);
-- point the next point to be processed.
newindex : integer := points'first-1;
pt : Link_to_Vector;
begin
if order
then pt := Head_Of(points(index));
else pt := Max_Extreme(points(index),n,-3,3);
Swap_to_Front(points(index),pt);
end if;
points(index) := Tail_Of(points(index));
for i in (index+1)..points'last loop
if not Is_Null(points(i))
then newindex := i;
end if;
exit when newindex >= points'first;
end loop;
if newindex < points'first
then for i in points'first..index loop
if not Is_Null(points(i))
then newindex := i;
end if;
exit when newindex >= points'first;
end loop;
end if;
index := newindex;
point := pt;
end Next_Point;
-- UPDATE ROUTINES :
procedure Merge ( mic : in out Mixed_Cell;
mixsub : in out Mixed_Subdivision ) is
tmp : Mixed_Subdivision := mixsub;
done : boolean := false;
begin
while not Is_Null(tmp) loop
declare
mic1 : Mixed_Cell := Head_Of(tmp);
last : List;
begin
if Equal(mic.nor,mic1.nor)
then for k in mic1.pts'range loop
last := mic1.pts(k);
while not Is_Null(Tail_Of(last)) loop
last := Tail_Of(last);
end loop;
Deep_Concat_Diff(mic.pts(k),mic1.pts(k),last);
end loop;
done := true;
else tmp := Tail_Of(tmp);
end if;
end;
exit when done;
end loop;
if done
then Deep_Clear(mic);
else Construct(mic,mixsub);
end if;
end Merge;
procedure Merge ( cells : in Mixed_Subdivision;
mixsub : in out Mixed_Subdivision ) is
tmp : Mixed_Subdivision := cells;
begin
while not Is_Null(tmp) loop
declare
mic : Mixed_Cell := Head_Of(tmp);
begin
Merge(mic,mixsub);
end;
tmp := Tail_Of(tmp);
end loop;
end Merge;
procedure Compute_New_Faces
( fs : in out Face_Structure; k,n : in natural;
point : in out Link_to_Vector; newfs : out Faces ) is
-- DESCRIPTION :
-- Given a point, the new faces will be computed and returned.
-- ON ENTRY :
-- fs a face structure;
-- k dimension of the faces to generate;
-- n dimension without the lifting;
-- point point that has to be in all the faces.
-- ON RETURN :
-- fs face structure with updated triangulation,
-- the new faces are not added, also the list is not updated;
-- point lifted conservatively w.r.t. fs.t;
-- newfs k-faces, which all contain the given point.
res,res_last : Faces;
procedure Append ( fc : in VecVec ) is
f : Face;
begin
f := new VecVec'(fc);
Append(res,res_last,f);
end Append;
procedure EnumLis is new Enumerate_Lower_Faces_in_List(Append);
procedure EnumTri is new Enumerate_Faces_in_Triangulation(Append);
begin
-- COMPUTE THE NEW FACES AND UPDATE fs :
if Is_Null(fs.t)
then
if Length_Of(fs.l) >= n
then
fs.t := Initial_Triangulation(n,fs.l,point);
if Is_Null(fs.t)
then EnumLis(fs.l,point,k);
else EnumTri(fs.t,point,k);
end if;
else
EnumLis(fs.l,point,k);
end if;
else
declare
newt : Triangulation;
begin
point(point'last) := Lift_to_Place(fs.t,point.all);
Update(fs.t,point,newt);
Enumtri(newt,point,k);
end;
end if;
Append(fs.l,fs.last,point);
newfs := res;
end Compute_New_Faces;
procedure Swap ( index : in natural; points : in out Array_of_Lists ) is
-- DESCRIPTION :
-- Swaps the first component of points with the component index.
tmp : List := points(index);
begin
points(index) := points(points'first);
points(points'first) := tmp;
end Swap;
procedure Swap ( index : in natural; mixsub : in out Mixed_Subdivision ) is
-- DESCRIPTION :
-- Swaps the first component of each cell with the component index.
tmp : Mixed_Subdivision := mixsub;
begin
while not Is_Null(tmp) loop
declare
mic : Mixed_Cell := Head_Of(tmp);
begin
Swap(index,mic.pts.all);
-- Set_Head(tmp,mic);
end;
tmp := Tail_Of(tmp);
end loop;
end Swap;
procedure Create_Cells
( index,n : in natural; mix : in Vector;
afs : in Array_of_Faces; lif : in Array_of_Lists;
nbsucc,nbfail : in out Standard_Floating_Vectors.Vector;
mixsub : out Mixed_Subdivision ) is
-- DESCRIPTION :
-- Creates a mixed subdivision by considering the faces
-- of component index first.
res : Mixed_Subdivision;
wrkmix : Vector(mix'range) := mix;
wrkafs : Array_of_Faces(afs'range) := afs;
wrklif : Array_of_Lists(lif'range) := lif;
begin
wrkmix(wrkmix'first) := mix(index); wrkmix(index) := mix(mix'first);
wrkafs(wrkafs'first) := afs(index); wrkafs(index) := afs(afs'first);
wrklif(wrklif'first) := lif(index); wrklif(index) := lif(lif'first);
Create_CS(n,wrkmix,wrkafs,wrklif,nbsucc,nbfail,res);
Swap(index,res);
mixsub := res;
end Create_Cells;
procedure Compute_New_Cells
( n : in natural; mix : in Vector; mic : in Mixed_Cell;
afs : in Array_of_Faces; index : in integer;
lifted : in Array_of_Lists; mixsub : out Mixed_Subdivision;
nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
-- DESCRIPTION :
-- Computes the new cells by considering a given mixed cell with
-- a tuple of faces.
-- ON ENTRY :
-- n dimension before lifting;
-- mix type of mixture;
-- mic a given mixed cell;
-- afs type of faces;
-- index component where new faces are computed for;
-- lifted tuple of lifted points;
-- nbsucc number of succesful tests of face-face combinations;
-- nbfail number of unsuccesful tests of face-face combinations.
-- ON RETURN :
-- mixsub the new mixed cells;
-- nbsucc updated number of succesful tests;
-- nbfail updated number of unsuccesful tests;
res,res2 : Mixed_Subdivision;
neiafs : Array_of_Faces(afs'range);
begin
neiafs(index) := Neighboring_Faces(mic,afs(index),index);
if not Is_Null(neiafs(index))
then
for i in neiafs'range loop
if i /= index
then neiafs(i) := Neighboring_Faces(mic,afs(i),i);
end if;
end loop;
if index /= 1
then Create_Cells(index,n,mix,neiafs,lifted,nbsucc,nbfail,res);
else Create_CS(n,mix,neiafs,lifted,nbsucc,nbfail,res);
end if;
for i in neiafs'range loop
if i /= index
then Shallow_Clear(neiafs(i));
end if;
end loop;
end if;
Shallow_Clear(neiafs(index));
res2 := Create(lifted,res); -- test on normals
Shallow_Clear(res);
mixsub := res2;
end Compute_New_Cells;
procedure Compute_New_Cells
( n : in natural; mix : in Vector;
mixsub : in Mixed_Subdivision; index : in integer;
newfaces : in Faces; fs : in Face_Structures;
newcells : out Mixed_Subdivision;
nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
res : Mixed_Subdivision;
afs : Array_of_Faces(fs'range);
lifted : Array_of_Lists(fs'range);
begin
-- CREATE NEW CELLS ONLY WITH THE NEW FACES :
for i in afs'range loop
if i = index
then afs(index) := newfaces;
else afs(i) := fs(i).f;
end if;
lifted(i) := fs(i).l;
end loop;
if not Is_Null(fs(index).t)
then
-- ENUMERATE ALL CELLS IN mixsub :
declare
tmp : Mixed_Subdivision := mixsub;
begin
while not Is_Null(tmp) loop
declare
newcells2 : Mixed_Subdivision;
begin
Compute_New_Cells(n,mix,Head_Of(tmp),afs,index,lifted,
newcells2,nbsucc,nbfail);
Merge(newcells2,res); Shallow_Clear(newcells2);
end;
tmp := Tail_Of(tmp);
end loop;
end;
else
-- COMPUTE ALL NEW CELLS AT ONCE :
declare
res2 : Mixed_Subdivision;
begin
if index /= 1
then Create_Cells(index,n,mix,afs,lifted,nbsucc,nbfail,res2);
else Create_CS(n,mix,afs,lifted,nbsucc,nbfail,res2);
end if;
res := Create(lifted,res2);
Shallow_Clear(res2);
end;
end if;
newcells := res;
end Compute_New_Cells;
-- BASIC VERSION : WITHOUT OUTPUT GENERICS :
procedure Dynamic_Lifting
( n : in natural; mix : in Vector;
points : in Array_of_Lists;
order,inter,conmv : in boolean; maxli : in natural;
mixsub : in out Mixed_Subdivision;
fs : in out Face_Structures;
nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
rest,inner : Array_of_Lists(points'range);
index,newindex : integer;
finished : boolean := false;
pt : Link_to_Vector;
nexli : natural := 1;
begin
Initialize(n,mix,points,mixsub,fs,rest,finished);
if not finished
then
index := First_Non_Empty(rest); newindex := index;
while not finished loop -- ITERATE UNTIL NO MORE POINTS :
Next_Point(n,rest,order,newindex,pt);
declare
point : Link_to_Vector := new vector(pt'first..pt'last+1);
newfa : Faces;
newcells : Mixed_Subdivision;
begin -- LIFT THE POINT CONSERVATIVELY :
point(pt'range) := pt.all;
point(point'last) := 1;
if inter and then Is_In(fs(index),point.all)
then
Clear(point); Construct(pt,inner(index));
else
nexli := Conservative_Lifting(mixsub,index,point.all);
if (maxli > 0) and then nexli > maxli
then Flatten(mixsub); Flatten(fs);
nexli := 1;
end if;
point(point'last) := nexli;
-- COMPUTE NEW FACES AND NEW CELLS :
Compute_New_Faces(fs(index),mix(index),n,point,newfa);
if not conmv or else not Zero_Contribution(n,fs,point.all,index)
then Compute_New_Cells(n,mix,mixsub,index,newfa,fs,
newcells,nbsucc,nbfail);
-- UPDATE THE MIXED SUBDIVISION AND THE FACE STRUCTURES :
Construct(newcells,mixsub); Shallow_Clear(newcells);
end if;
Construct(fs(index).f,newfa); Shallow_Clear(newfa);
end if;
end;
index := newindex;
finished := (index < rest'first);
end loop;
end if;
if inter -- lift out the interior points
then for i in inner'range loop
Constant_Lifting(inner(i),nexli,fs(i).l,fs(i).last);
Shallow_Clear(inner(i));
end loop;
end if;
exception
when numeric_error -- probably due to a too high lifting
=> Flatten(mixsub); Flatten(fs);
Dynamic_Lifting(n,mix,rest,order,inter,conmv,maxli,mixsub,fs,
nbsucc,nbfail);
end Dynamic_Lifting;
-- EXTENDED VERSIONS : WITH OUTPUT GENERICS :
procedure Dynamic_Lifting_with_Flat
( n : in natural; mix : in Vector; points : in Array_of_Lists;
order,inter,conmv : in boolean; maxli : in natural;
mixsub : in out Mixed_Subdivision;
fs : in out Face_Structures;
nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
rest,inner : Array_of_Lists(points'range);
index,newindex : integer;
finished : boolean := false;
pt : Link_to_Vector;
nexli : natural := 1;
begin
Initialize(n,mix,points,mixsub,fs,rest,finished);
if not finished
then
index := First_Non_Empty(rest); newindex := index;
while not finished loop -- ITERATE UNTIL NO MORE POINTS :
Next_Point(n,rest,order,newindex,pt);
declare
point : Link_to_Vector := new vector(pt'first..pt'last+1);
newfa : Faces;
newcells : Mixed_Subdivision;
begin -- LIFT THE POINT CONSERVATIVELY :
point(pt'range) := pt.all;
point(point'last) := 1;
if inter and then Is_In(fs(index),point.all)
then
Clear(point); Construct(pt,inner(index));
else
nexli := Conservative_Lifting(mixsub,index,point.all);
if (maxli > 0) and then nexli > maxli
then Before_Flattening(mixsub,fs);
Flatten(mixsub); Flatten(fs);
nexli := 1;
end if;
point(point'last) := nexli;
-- COMPUTE NEW FACES AND NEW CELLS :
Compute_New_Faces(fs(index),mix(index),n,point,newfa);
if not conmv or else not Zero_Contribution(n,fs,point.all,index)
then Compute_New_Cells(n,mix,mixsub,index,newfa,fs,
newcells,nbsucc,nbfail);
-- UPDATE THE MIXED SUBDIVISION AND THE FACE STRUCTURES :
Construct(newcells,mixsub); Shallow_Clear(newcells);
end if;
Construct(fs(index).f,newfa); Shallow_Clear(newfa);
end if;
end;
index := newindex;
finished := (index < rest'first);
end loop;
end if;
if inter -- lift out the interior points
then for i in inner'range loop
Constant_Lifting(inner(i),nexli,fs(i).l,fs(i).last);
Shallow_Clear(inner(i));
end loop;
end if;
exception
when numeric_error -- probably due to a too high lifting
=> Before_Flattening(mixsub,fs);
Flatten(mixsub); Flatten(fs);
Dynamic_Lifting_with_Flat(n,mix,rest,order,inter,conmv,maxli,mixsub,
fs,nbsucc,nbfail);
end Dynamic_Lifting_with_Flat;
procedure Dynamic_Lifting_with_New
( n : in natural; mix : in Vector; points : in Array_of_Lists;
order,inter,conmv : in boolean; maxli : in natural;
mixsub : in out Mixed_Subdivision;
fs : in out Face_Structures;
nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
rest,inner : Array_of_Lists(points'range);
index,newindex : integer;
finished : boolean := false;
pt : Link_to_Vector;
nexli : natural := 1;
begin
Initialize(n,mix,points,mixsub,fs,rest,finished);
Process_New_Cells(mixsub,0,Head_Of(fs(fs'last).last).all);
if not finished
then
index := First_Non_Empty(rest); newindex := index;
while not finished loop -- ITERATE UNTIL NO MORE POINTS :
Next_Point(n,rest,order,newindex,pt);
declare
point : Link_to_Vector := new vector(pt'first..pt'last+1);
newfa : Faces;
newcells : Mixed_Subdivision;
begin -- LIFT THE POINT CONSERVATIVELY :
point(pt'range) := pt.all;
point(point'last) := 1;
if inter and then Is_In(fs(index),point.all)
then
Clear(point); Construct(pt,inner(index));
else
nexli := Conservative_Lifting(mixsub,index,point.all);
if (maxli > 0) and then nexli > maxli
then Flatten(mixsub); Flatten(fs);
nexli := 1;
end if;
point(point'last) := nexli;
-- COMPUTE NEW FACES AND NEW CELLS :
Compute_New_Faces(fs(index),mix(index),n,point,newfa);
if not conmv or else not Zero_Contribution(n,fs,point.all,index)
then Compute_New_Cells(n,mix,mixsub,index,newfa,fs,newcells,
nbsucc,nbfail);
Process_New_Cells(newcells,index,point.all);
-- UPDATE THE MIXED SUBDIVISION AND THE FACE STRUCTURES :
Construct(newcells,mixsub); Shallow_Clear(newcells);
end if;
Construct(fs(index).f,newfa); Shallow_Clear(newfa);
end if;
end;
index := newindex;
finished := (index < rest'first);
end loop;
end if;
if inter -- lift out the interior points
then for i in inner'range loop
Constant_Lifting(inner(i),nexli,fs(i).l,fs(i).last);
Shallow_Clear(inner(i));
end loop;
end if;
exception
when numeric_error -- probably due to a too high lifting
=> Flatten(mixsub); Flatten(fs);
Dynamic_Lifting_with_New(n,mix,rest,order,inter,conmv,maxli,mixsub,fs,
nbsucc,nbfail);
end Dynamic_Lifting_with_New;
procedure Dynamic_Lifting_with_Flat_and_New
( n : in natural; mix : in Vector; points : in Array_of_Lists;
order,inter,conmv : in boolean; maxli : in natural;
mixsub : in out Mixed_Subdivision;
fs : in out Face_Structures;
nbsucc,nbfail : in out Standard_Floating_Vectors.Vector ) is
rest,inner : Array_of_Lists(points'range);
index,newindex : integer;
finished : boolean := false;
pt : Link_to_Vector;
nexli : natural := 1;
begin
Initialize(n,mix,points,mixsub,fs,rest,finished);
Process_New_Cells(mixsub,0,Head_Of(fs(fs'last).last).all);
if not finished
then
index := First_Non_Empty(rest); newindex := index;
while not finished loop -- ITERATE UNTIL NO MORE POINTS
Next_Point(n,rest,order,newindex,pt);
declare
point : Link_to_Vector := new vector(pt'first..pt'last+1);
newfa : Faces;
newcells : Mixed_Subdivision;
begin -- LIFT THE POINT CONSERVATIVELY :
point(pt'range) := pt.all;
point(point'last) := 1;
if inter and then Is_In(fs(index),point.all)
then
Clear(point); Construct(pt,inner(index));
else
nexli := Conservative_Lifting(mixsub,index,point.all);
if (maxli > 0) and then nexli > maxli
then Before_Flattening(mixsub,fs);
Flatten(mixsub); Flatten(fs);
nexli := 1;
end if;
point(point'last) := nexli;
-- COMPUTE NEW FACES AND NEW CELLS :
Compute_New_Faces(fs(index),mix(index),n,point,newfa);
if not conmv or else not Zero_Contribution(n,fs,point.all,index)
then Compute_New_Cells(n,mix,mixsub,index,newfa,fs,newcells,
nbsucc,nbfail);
Process_New_Cells(newcells,index,point.all);
-- UPDATE THE MIXED SUBDIVISION AND THE FACE STRUCTURES :
Construct(newcells,mixsub); Shallow_Clear(newcells);
end if;
Construct(fs(index).f,newfa); Shallow_Clear(newfa);
end if;
end;
index := newindex;
finished := (index < rest'first);
end loop;
end if;
if inter -- lift out the interior points
then for i in inner'range loop
Constant_Lifting(inner(i),nexli,fs(i).l,fs(i).last);
Shallow_Clear(inner(i));
end loop;
end if;
exception
when numeric_error -- probably due to a too high lifting
=> Before_Flattening(mixsub,fs);
Flatten(mixsub); Flatten(fs);
Dynamic_Lifting_with_Flat_and_New
(n,mix,rest,order,inter,conmv,maxli,mixsub,fs,nbsucc,nbfail);
end Dynamic_Lifting_with_Flat_and_New;
end Dynamic_Mixed_Subdivisions;