Artifact 44f1f4147dd9e1cf35cc69fdcc877e8dd75abac93ab0115dc8ad8f773084d247:
- Executable file
r37/packages/plot/plotimp2.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 15090) [annotate] [blame] [check-ins using] [more...]
module plotimp2; % Implicit plot: compute the varity {x,y|f(x,y)=c}. % Author: Herbert Melenk, ZIB Berlin. % data structure: % % point = {x,y,f(x,y),t1,t2,t3,...} % where ti are the numbers of the triangles which % have this point as a node. % The point list is unique - adjacent triangles share % the list for equal nodes. The node numbers are % updated in place. % % triangle = {t#,p1,p2,p3} % triangles are stored in triangle vector by number % fluid '(imp2!-triacount!* imp2!-trias!* !*imp2!-trace); fluid '(!*show_grid !*test_plot plot!-contour* plot!-refine!*); imp2!-triacount!*:=0; symbolic procedure ploteval2xyimpl(rx,ry,f,x,y); begin scalar ll,l,form,g; for each c in plot!-contour!* do << form := plot!-form!-prep({'difference,cdr f,c},x,y); l:=imp2!-plot(car rx,cadr rx, car ry,cadr ry, plot!-points(nil),form); ll:=append(ll,l); >>; if !*show_grid and null cdr plot!-contour!* then g:= imp2!-show!-meshes(); plotdriver(plot!-2imp,x,y,ll,g,car rx,cadr rx,car ry,cadr ry); end; symbolic procedure imp2!-init(); << imp2!-finit(); if null imp2!-trias!* then imp2!-trias!* :=mkxvect()>>; symbolic procedure imp2!-finit(); <<if imp2!-trias!* then for i:=0:imp2!-triacount!* do xputv(imp2!-trias!*,i,nil); imp2!-triacount!*:=0; >>; symbolic procedure mk!-point(x0,y0,fcn); {x0,y0,apply2(fcn,x0,y0)}; !#if (member 'csl lispsystem!*) symbolic procedure deletip1 (u,v); % Auxiliary function for DeletIP. pairp cdr v and (if u=cadr v then rplacd(v,cddr v) else deletip1(u,cdr v)); symbolic procedure deletip (u,v); % Destructive DELETE. if not pairp v then v else if u=car v then cdr v else <<deletip1(u,v); v>>; !#endif symbolic procedure imp2!-delete!-pt!-reference(i,p); cdr cddr p := deletip(i,cdddr p); symbolic procedure mk!-tria(i,p1,p2,p3); % make a triangle from 3 points. If the number is given, % reuse it. Otherwise generate a fresh number. begin scalar p; integer i; i := i or (imp2!-triacount!* := imp2!-triacount!* #+1); p:={i,p1,p2,p3,imp2!-tria!-zerop!*(p1,p2,p3)}; xputv(imp2!-trias!*,i,p); aconc(p1,i); aconc(p2,i); aconc(p3,i); if !*imp2!-trace then print!-tria("new triangle ",p); return p; end; symbolic procedure print!-tria(tx,w); <<prin2 tx; prin2 car w; w:=cdr w; prin2l {{car car w,cadr car w,{caddr car w}}, {car cadr w,cadr cadr w,{caddr cadr w}}, {car caddr w,cadr caddr w,{caddr caddr w}}}; terpri(); >>; symbolic procedure imp2!-tria!-zerop!*(p1,p2,p3); % Here I test whether the triangle contains a zero line. begin scalar f1,f2,f3; f1 := caddr p1; f2 := caddr p2; f3 := caddr p3; return f1*f2 <= 0.0 or f1*f3 <= 0.0; end; symbolic procedure imp2!-tria!-zerop(w); % Fast access to stored zerop property. cadddr cdr w; symbolic procedure imp2!-neighbours p; % Compute the direct neighbours of p - the traingles which have % two nodes respectively one complete edge in common with p. begin scalar l,p1,p2,p3; integer n; if fixp p then p:=xgetv(imp2!-trias!*,p); n:= car p; p:=cdr p; p1:=delete(n,cdddr car p); p2:=delete(n,cdddr cadr p); p3:=delete(n,cdddr caddr p); l:={find!-one!-common(p1,p2), find!-one!-common(p2,p3), find!-one!-common(p3,p1)}; while nil memq l do l:=deletip(nil,l); return for each w in l collect xgetv(imp2!-trias!*,w); end; symbolic procedure imp2!-edge!-length(p1,p2); begin scalar dx,dy; fdeclare('imp2!-edge!-length,dx,dy); dx:=thefloat car p1 - thefloat car p2; dy:=thefloat cadr p1 - thefloat cadr p2; return sqrt(dx*dx + dy*dy) end; symbolic procedure imp2!-tria!-surface w; begin scalar x1,x2,x3,y1,y2,y3,p1,p2,p3; fdeclare('imp2!-tria!-surface,x1,x2,x3,y1,y2,y3); w:=cdr w; x1:=car (p1:=car w); y1:=cadr p1; x2:=car (p2:=cadr w); y2:=cadr p2; x3:=car (p3:=caddr w); y3:=cadr p3; return abs ((0.5*(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2)))); end; symbolic procedure imp2!-tria!-length w; begin scalar p1,p2,p3; w:=cdr w; p1:=car w; p2:=cadr w; p3:=caddr w; return (0.3*(imp2!-edge!-length(p1,p2) + imp2!-edge!-length(p2,p3) + imp2!-edge!-length(p3,p1))); end; symbolic procedure imp2!-tria!-midpoint(w); <<w:= cdr w; {(thefloat car car w + thefloat car cadr w + thefloat car caddr w)/3.0, (thefloat cadr car w + thefloat cadr cadr w + thefloat cadr caddr w)/3.0} >>; symbolic procedure imp2!-tria!-goodpoint(w,fn); % Here I assume that there is a zero in the triangle. Compute % a point inside the triangle which is as close as possible % to the desired line, but without local recomputation of % function values. begin scalar p1,p2,p3,f1,f2,f3,s1,s2,s3; w:=cdr w; p1:=car w; p2:=cadr w; p3:=caddr w; if (f1:=caddr p1)=0.0 then return {car p1,cadr p1} else if (f2:=caddr p2)=0.0 then return {car p2,cadr p2} else if (f3:=caddr p3)=0.0 then return {car p3,cadr p3}; s1:=f1<0.0; s2:=f2<0.0; s3:=f3<0.0; return if s1=s2 then imp2!-tria!-goodpoint1(p1,f1,p3,f3,p2,f2,fn) else if s1=s3 then imp2!-tria!-goodpoint1(p1,f1,p2,f2,p3,f3,fn) else imp2!-tria!-goodpoint1(p2,f2,p1,f1,p3,f3,fn) end; %symbolic procedure imp2!-tria!-goodpoint1(p1,f1,p2,f2,p3,f3,fn); % % Now I know that f2 has the opposite sign to f1 and f3. % % I take the linearly interpolated zero of f on p1-p2 and p2-p3 % % return their arithmetic mean value which lies inside the % % triangle. % begin scalar x1,x2,y1,y2,s; % fdeclare (x1,x2,y1,y2,s,f1,f2,f3); % s:=f1-f2; % x1:=(f1*thefloat car p2 - f2*thefloat car p1)/s; % y1:=(f1*thefloat cadr p2 - f2*thefloat cadr p1)/s; % s:=f3-f2; % x2:=(f3*thefloat car p2 - f2*thefloat car p3)/s; % y2:=(f3*thefloat cadr p2 - f2*thefloat cadr p3)/s; % return {(x1+x2)*0.5, (y1+y2)*0.5}; % end; symbolic procedure imp2!-tria!-goodpoint1(p1,f1,p2,f2,p3,f3,fn); % Now I know that f2 has the opposite sign to f1 and f3. % F1 and f3 are non-zero. % I use the midpoint of the p1-p3 edge and look for an % approximation of a zero on the connection of the midpoint % and p2 using regula falsi. begin scalar x1,x2,y1,y2,x3,y3,s; fdeclare (x1,x2,x3,y1,y2,y3,s,f1,f2,f3); f1:=(f1+f3)*0.5; x1:=(thefloat car p1 + thefloat car p3)*0.5; y1:=(thefloat cadr p1 + thefloat cadr p3)*0.5; x2:=car p2; y2:=cadr p2; s:=f2-f1; x3:=(x1*f2 - x2*f1)/s; y3:=(y1*f2 - y2*f1)/s; f3:=apply2(fn,x3,y3); if f2*f3>=0 then <<s:=f1-f3; x3:=(x3*f1-x1*f3)/s; y3:=(y3*f1-y1*f3)/s>> else <<s:=f2-f3; x3:=(x3*f2-x2*f3)/s; y3:=(y3*f2-y2*f3)/s>>; done: return{x3,y3}; end; symbolic procedure imp2!-tri!-refine!-one!-tria (w,fn); % Refine one triangle by inserting a new point in the mid % of the longest edge. Also, refine the triangle adjacent % to that edge with the same point. begin scalar p1,p2,p3,pn,xn,yn,new,nb,y; integer i; scalar x1,x2,y1,y2,d1,d2,d3,s; fdeclare (x1,x2,y1,y2,s,d1,d2,d3); if fixp w then w :=xgetv(imp2!-trias!*,w); % record. if !*imp2!-trace then print!-tria("refine ",w); i:=car w; w :=cdr w; % delete reference to this triangle. imp2!-delete!-pt!-reference(i,car w); imp2!-delete!-pt!-reference(i,cadr w); imp2!-delete!-pt!-reference(i,caddr w); % find longest edge d1:=imp2!-edge!-length(car w,cadr w); d2:=imp2!-edge!-length(cadr w,caddr w); d3:=imp2!-edge!-length(car w,caddr w); if d1>=d2 and d1>=d3 then <<p1:=car w; p2:=cadr w; p3:=caddr w>> else if d2>=d1 and d2>=d3 then <<p1:=cadr w; p2:=caddr w; p3:=car w>> else <<p1:=car w; p2:=caddr w, p3:=cadr w>>; % identify the neighbour triangle. nb := find!-one!-common(cdddr p1,cdddr p2); % compute new point almost at the mid. s:=0.45+0.01*random(10); x1:=car p1; y1:=cadr p1; x2:=car p2; y2:=cadr p2; xn:=x1*s+x2*(1.0-s); yn:=y1*s+y2*(1.0-s); pn:=mk!-point(xn,yn,fn); construct: % construct new triangles new:=mk!-tria(i,p1,pn,p3).new; new:=mk!-tria(nil,p2,pn,p3).new; % update the neighbour, if there is one if null nb then return new; w:=nb; nb:=nil; if fixp w then w :=xgetv(imp2!-trias!*,w); % record. i:=car w; w:=cdr w; imp2!-delete!-pt!-reference(i,car w); imp2!-delete!-pt!-reference(i,cadr w); imp2!-delete!-pt!-reference(i,caddr w); % find the far point p3 := if not((y:=car w) eq p1 or y eq p2) then y else if not((y:=cadr w) eq p1 or y eq p2) then y else caddr w; goto construct; end; %symbolic procedure imp2!-tri!-refine!-one!-tria!-center (w,fn); % % Refine one triangle by inserting a new point in the center. % begin scalar p1,p2,p3,pn,xn,yn,new,nb,y; integer i; % scalar x1,x2,x3,y1,y2,y3; % fdeclare (x1,x2,x3,y1,y2,y3); % if fixp w then w :=xgetv(imp2!-trias!*,w); % record. % if !*imp2!-trace then print!-tria("refine ",w); % i:=car w; w :=cdr w; % % % delete reference to this triangle. % imp2!-delete!-pt!-reference(i,car w); % imp2!-delete!-pt!-reference(i,cadr w); % imp2!-delete!-pt!-reference(i,caddr w); % % % compute center. % p1:=car w; p2:=cadr w; p3:=caddr w; % x1:=car p1; y1:=cadr p1; % x2:=car p2; y2:=cadr p2; % x3:=car p3; y3:=cadr p3; % xn:=(x1+x2+x3)*0.33333; % yn:=(y1+y2+y3)*0.33333; % pn:=mk!-point(xn,yn,fn); %construct: % % construct new triangles % new:=mk!-tria(i,p1,p2,pn).new; % new:=mk!-tria(nil,p2,p3,pn).new; % new:=mk!-tria(nil,p1,p3,pn).new; % return new; % end; symbolic procedure find!-one!-common(u,v); % fast equivalent to "car intersection(u,v)". if null u then nil else if memq(car u,v) then car u else find!-one!-common(cdr u,v); %%%%%% Main program for implicit plot. symbolic procedure imp2!-plot(xlo,xhi,ylo,yhi,pts,fcn); begin scalar p1,p2,p3,p4,sf,z; imp2!-init(); % setup initial triangularization. p1:=mk!-point(xlo,ylo,fcn); p2:=mk!-point(xhi,ylo,fcn); p3:=mk!-point(xhi,yhi,fcn); p4:=mk!-point(xlo,yhi,fcn); mk!-tria(nil,p1,p2,p3); mk!-tria(nil,p1,p3,p4); sf:=((xhi-xlo)+(yhi-ylo))*1.5/float pts; z:=imp2!-plot!-refine(sf,fcn); if !*imp2!-trace then for each w in z do print!-tria("zero triangle:",w); if !*test_plot then print "collect"; z:=imp2!-plot!-collect(z); if !*test_plot then print "draw"; z:=imp2!-plot!-draw(z,fcn); if not !*show_grid then imp2!-finit(); return z; end; symbolic procedure imp2!-plot!-refine(sflimit,fcn); begin scalar z,zn,w,s,limit2,again; integer i,j; limit2 := 0.15*sflimit; % In the first stage I refine all triangles until they % are fine enough for a coarse overview. again := t; if !*test_plot then print "phase1"; phase1: z:=nil; again:=nil; for i:=imp2!-triacount!* step -1 until 1 do << w := xgetv(imp2!-trias!*,i); if imp2!-tria!-length w > sflimit then <<again:=t; imp2!-tri!-refine!-one!-tria (w,fcn)>> else if not again and imp2!-tria!-zerop w then z:=car w.z; >>; j:=j #+ 1; if j+j #< plot!-refine!* or again then goto phase1; % Next I refine further only the triangles which contain a zero. % I must store in z only the numbers of triangles because these % may be modified as side effects by copying. if !*test_plot then print "phase2"; phase2: for each w in z do if (s:=imp2!-tria!-length (w:=xgetv(imp2!-trias!*,w))) >limit2 then <<for each q in imp2!-tri!-refine!-one!-tria (w,fcn) do if imp2!-tria!-zerop q and not memq(car q,zn) then zn:=car q.zn>>; z:=zn; zn:=nil; if z then goto phase2; % In the third phase I refine those critical triangles futher: % non-zeros with two zero neighbours. These may be turning points % or bifurcations. if !*test_plot then print "phase3"; phase3: for i:=imp2!-triacount!* step -1 until 1 do imp2!-refine!-phase3(i,i,plot!-refine!*/2,fcn,limit2*0.3); % Return the final list of zeros in ascending order. z:=for i:=1:imp2!-triacount!* join if imp2!-tria!-zerop(w:=xgetv(imp2!-trias!*,i)) then {w}; return z; end; symbolic procedure imp2!-refine!-phase3(i,low,lev,fn,lth); begin scalar w; integer c; if lev=0 then return nil; w:=if numberp i then xgetv(imp2!-trias!*,i) else i; if car w<low or imp2!-tria!-length w < lth then return nil; lev:=lev-1; for each q in imp2!-neighbours w do if imp2!-tria!-zerop q then c:=c+1; if imp2!-tria!-zerop w then (if c eq 2 then return nil) else (if c #< 2 then return nil); for each q in imp2!-tri!-refine!-one!-tria (w,fn) do imp2!-refine!-phase3(q,low,lev,fn,lth); end; symbolic procedure imp2!-plot!-collect(z); begin scalar lines,l,q,p,s; for each w1 in z do for each w2 in imp2!-neighbours w1 do if car w2 > car w1 and imp2!-tria!-zerop w2 then q:=(w1.w2) . q; lineloop: if null q then return lines; l:={caar q, (p:=cdar q)}; q:= cdr q; % first I extend the back end. while q and p do << if(s:= atsoc(p,q)) then l:=nconc(l,{p:=cdr s}) else if(s:=rassoc(p,q)) then l:=nconc(l,{p:=car s}); if s then q:=delete(s,q) else p:=nil; >>; % next I extend the front end. p:=car l; while q and p do << if(s:=rassoc(p,q)) then l:=(p:=car s).l else if(s:= atsoc(p,q)) then l:=(p:=cdr s).l; if s then q:=delete(s,q) else p:=nil; >>; lines := l.lines; goto lineloop; end; symbolic procedure imp2!-plot!-draw(l,fn); begin scalar r,s,q; q:=for each w in l collect <<r:=nil; for each q in w join <<s:=imp2!-tria!-goodpoint(q,fn); if r neq s then {r:=s}>> >>; return q; end; symbolic procedure imp2!-show!-meshes(); % generate plot of meshes; begin scalar d,l,w,s,p1,p2; integer i; i:=1; loop: w:=xgetv(imp2!-trias!*,i); if null w then <<imp2!-finit(); return l>>; w:=cdr w; d:= {{car car w, cadr car w}, {car cadr w,cadr cadr w}, {car caddr w,cadr caddr w}, {car car w, cadr car w}} ; while d and cdr d do <<p1:=car d; p2:=cadr d; d:=cdr d; if car p1 > car p2 then <<w:=p2;p2:=p1;p1:=w>>; s:={p1,p2}; if not member(s,l) then l:=s.l >>; i:=i+1; goto loop; end; endmodule; % plotimp2 end;