Unit Grafisch; INTERFACE Uses Omcirkel, Unit2, Graphics; type row4double = array[1..4] of double; procedure HdBcheck(naam : string); procedure HdBpost1(naam,kenmerk : string ; lengte : integer); procedure HdBpost2(kenmerk : string ; nivos : integer); procedure HdBpost3(nivos : integer); procedure HdBpost4(lengte : integer); IMPLEMENTATION type rijtje = array[1..3] of byte; {************************************} procedure permutatie(L, n : byte ; var r : rijtje); { Permutations of array r(n) -------------------------- L := ordinal number of permutation (input) n := number of elements in: (input) r := row containing the permutaton (output) } var h : integer; k, p, q, s, m : byte; begin h := 1 ; for k := 1 to n do begin r[k] := k ; h := h*k end; { Which item: } for k := 1 to n do begin p := n-k+1 ; h := h div p ; q := (((L-1) div h) mod p) + 1 ; if q = 1 then Continue; { Otherwise rotate data: } s := r[k+q-1] ; for m := (k+q-1) downto (k+1) do r[m] := r[m-1] ; r[k] := s ; end; end; {************************************} function log2(n : integer) : byte; { Logarithm base 2 ---------------- } var k, L : byte; begin log2 := 0; if n <= 0 then Exit; for k := 1 to 32 do begin n := n div 2 ; L := k; if n = 0 then Break ; end; log2 := L; end; {************************************} procedure Contour_Vierhoek(meest : byte ; xq, yq, zq : row4double); var i, j, k, n, nq1, nq2, nop : byte; xd, yd, zd : array[1..3] of double ; qmin, qmax, zero : double; {************************************} procedure onder_boven(q : double ; most : byte ; var n1, n2 : byte); { Determine lower & upper boundary values --------------------------------------- } var m, grens, k : byte; begin n1 := 0 ; n2 := most ; m := most ; grens := log2(most) ; for k := 1 to grens do begin m := m div 2 ; n1 := n1 + m ; n2 := n2 - m ; if q <= n1/most then n1 := n1 - m ; if q >= n2/most then n2 := n2 + m ; end; end; {************************************} procedure Contour_Driehoek; { Contourlines in a Triangle -------------------------- } const eerst : boolean = true; np : array[1..3,1..6] of byte = ((0,0,0,0,0,0),(0,0,0,0,0,0),(0,0,0,0,0,0)); klas : array[1..27] of byte = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0); move : array[1..27] of byte = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0); var Hoog : integer; fd : array[1..3] of double; k, more, merk, keus : byte; np1, np2, np3 : byte; x1, x2, x3 : double; y1, y2, y3 : double; f1, f2, f3 : double; noemer, xc, yc : double ; {************************************} procedure Setup_Contours; { Initialize np[3,6],klas[27],move[27] ------------------------------------ All possible cases to be classified: === ==< ==> =<= =<< =<> =>= =>< =>> <== <=< <=> <<= <<< <<> <>= <>< <>> >== >=< >=> ><= ><< ><> >>= >>< >>> Being equivalent to counting base 3, where 0 : = , 1 : < , 2 : > . Resulting in 6 equivalence classes: 1 2 3 4 5 6 keer === ==< =<< =<> <<< <<> 1 ==> =>> =>< >>> >>< 2 =<= <=< <=> <>< 1 =>= >=> >=< ><> 2 <== <<= <>= ><< 1 >== >>= ><= <>> 2 } var L, i, j, n, m, keer : byte; num, rij : rijtje; begin for L := 1 to 6 do begin permutatie(L, 3, rij) ; for j := 1 to 3 do begin np[j,L] := rij[j] ; end; end; for L := 1 to 27 do begin klas[L] := 0 ; move[L] := 0 ; end; { Equivalence classes: } n := 0 ; for L := 1 to 27 do begin if klas[L] > 0 then Continue; n := n+1 ; m := L-1 ; { Decimals base 3: } for i := 1 to 3 do begin num[3-i+1] := m mod 3 ; m := m div 3 ; end; for keer := 1 to 2 do begin for i := 1 to 6 do begin { Inverse permutation: } for j := 1 to 3 do rij[np[j,i]] := j ; m := 0 ; for j := 1 to 3 do m := m*3+num[rij[j]] ; m := m+1 ; { Assign equivalent elements: } if klas[m] = 0 then klas[m] := n ; { Class } if move[m] = 0 then move[m] := i ; { Permutation } end; for i := 1 to 3 do num[i] := (3-num[i]) mod 3 ; m := 0 ; for j := 1 to 3 do m := m*3 + num[rij[j]] ; m := m+1 ; if klas[m] > 0 then Break ; end; { < and > exchanged 2nd time } end; end; {************************************} begin if eerst then begin Setup_Contours ; eerst := false; end; Hoog := Form2.Image1.Height; for k := 1 to 3 do fd[k] := zd[k] - zero ; { Coding each case: } more := 0 ; for k := 1 to 3 do begin if fd[k] = 0 then more := more*3 ; if fd[k] < 0 then more := more*3+1 ; if fd[k] > 0 then more := more*3+2 ; end; merk := klas[more+1] ; if (merk = 1) or (merk = 3) or (merk = 5) then Exit ; { === =<< >>> } { Permutation of (1,2,3): } keus := move[more+1] ; np1 := np[1,keus] ; np2 := np[2,keus] ; np3 := np[3,keus] ; x1 := xd[np1] ; x2 := xd[np2] ; x3 := xd[np3] ; y1 := yd[np1] ; y2 := yd[np2] ; y3 := yd[np3] ; f1 := fd[np1] ; f2 := fd[np2] ; f3 := fd[np3] ; Case merk of 1, 3, 5 : Exit ; { === , =<< , <<< } 2 : { line through side: ==< } with Form2.Image1.Canvas do begin MoveTo(Trunc(x1),Trunc(Hoog-y1)) ; LineTo(Trunc(x2),Trunc(Hoog-y2)) ; Exit; end; 4 : { line through corner: =<> } with Form2.Image1.Canvas do begin MoveTo(Trunc(x1),Trunc(Hoog-y1)) ; noemer := 1/(f2-f3) ; xc := (f2*x3-f3*x2)*noemer; yc := (f2*y3-f3*y2)*noemer; LineTo(Trunc(xc),Trunc(Hoog-yc)) ; Exit; end; 6 : { line through two sides: <<> } with Form2.Image1.Canvas do begin noemer := 1/(f1-f3) ; xc := (f1*x3-f3*x1)*noemer; yc := (f1*y3-f3*y1)*noemer; MoveTo(Trunc(xc),Trunc(Hoog-yc)) ; noemer := 1/(f2-f3) ; xc := (f2*x3-f3*x2)*noemer; yc := (f2*y3-f3*y2)*noemer; LineTo(Trunc(xc),Trunc(Hoog-yc)) ; Exit; end; end; end; {************************************} begin qmin := 1 ; qmax := 0 ; for i := 1 to 4 do begin if zq[i] < qmin then qmin := zq[i] ; if zq[i] > qmax then qmax := zq[i] ; end; onder_boven(qmin, meest, nq1, nop) ; onder_boven(qmax, meest, nop, nq2) ; { Loop through Triangles: } xd[3] := 0.25*(xq[1]+xq[2]+xq[3]+xq[4]); yd[3] := 0.25*(yq[1]+yq[2]+yq[3]+yq[4]); zd[3] := 0.25*(zq[1]+zq[2]+zq[3]+zq[4]); for i := 1 to 4 do begin for j := 1 to 2 do begin k := ((i+j-2) mod 4) + 1; xd[j] := xq[k] ; yd[j] := yq[k] ; zd[j] := zq[k] ; end; { Loop through Levels: } for n := nq1 to nq2 do begin zero := n/meest ; { Elementary contours: } Contour_Driehoek; end; end; end; {************************************} procedure HdBcontour(most : byte ; naam : string); var i : byte; file1 : text; MaxY, limc, L, grens, nel : integer; Xbot, Xtop, Ybot, Ytop, Fbot, Ftop : double; xq, yq, zq : row4double; ff : double; b : array of double; begin limc := ndr*nph; SetLength(b,limc+1); Assign(file1, naam); Reset(file1); L := 0; while not EoF(file1) do begin L := L + 1; Readln(file1, ff); b[L] := ff; end; Close(file1); limc := L; Xbot := xx(1) ; Xtop := xx(1) ; for L := 2 to limc do begin if Xtop < xx(L) then Xtop := xx(L) ; if Xbot > xx(L) then Xbot := xx(L) ; end; Ybot := yy(1) ; Ytop := yy(1) ; for L := 2 to limc do begin if Ytop < yy(L) then Ytop := yy(L) ; if Ybot > yy(L) then Ybot := yy(L) ; end; Fbot := b[1] ; Ftop := b[1] ; for L := 2 to limc do begin if Ftop < b[L] then Ftop := b[L] ; if Fbot > b[L] then Fbot := b[L] ; end; { Loop through mesh: } MaxY := Form2.Image1.Height; grens := (ndr-1)*(nph-1); for nel := 1 to grens do begin PoolWeb(nel, kon[1], kon[2], kon[3], kon[4]); { At quadrilaterals, normed: } for i := 1 to 4 do begin L := kon[i] ; xq[i] := MaxY*(xx(L)-Xbot)/(Xtop-Xbot) ; yq[i] := MaxY*(yy(L)-Ybot)/(Ytop-Ybot) ; zq[i] := (b[L]-Fbot)/(Ftop-Fbot) ; end; Contour_Vierhoek(most, xq, yq, zq); end; SetLength(b,0); end; {************************************} procedure VectorPlot(naam : string; lengte : integer); { Vectorplot velocities --------------------- } var Breed : integer; x, y : integer; u, v : double; k, n1, n2 : integer; file1 : text; begin Breed := Form2.Image1.Height; Form2.Image1.Canvas.Pen.Color := clGreen ; Assign(file1, naam); Reset(file1); k := 0; while not EoF(file1) do begin Readln(file1, u); Readln(file1, v); k := k + 1; Midden_Knopen(k, n1,n2); x := Trunc(Breed*0.5*0.1*(xx(n1)+xx(n2))); y := Trunc(Breed*(1-0.5*0.1*(yy(n1)+yy(n2)))); Form2.Image1.Canvas.MoveTo(x,y); x := x + Trunc(Breed*0.01*lengte*u); y := y - Trunc(Breed*0.01*lengte*v); Form2.Image1.Canvas.LineTo(x,y); end; Close(file1); end; {************************************} procedure Cirkel_Analytisch(nul : vier_ints); { Plot exact velocities on circle ------------------------------- } var phi, uu, vv, faktor : double; xo, yo, sx, sy : integer; xn, yn : integer; i, limiet : byte; begin xo := nul.xo ; yo := nul.yo ; sx := nul.sx ; sy := nul.sy ; xn := xo+sx ; yn := yo-sy ; Form2.Image1.Canvas.MoveTo(xn,yn) ; Form2.Image1.Canvas.Pen.Color := clBlue; ; limiet := 25 ; { Horizontal component ... } faktor := 0.5*Pi/(limiet-1) ; for i := 1 to limiet do begin phi := faktor*(i-1) ; uu := 1-cos(2*phi) ; xn := xo+Trunc(cos(phi)*sx) ; yn := yo-Trunc(sy*(1+uu)) ; Form2.Image1.Canvas.LineTo(xn,yn) ; end; { Vertical component ... } xn := xo+sx ; yn := yo-sy ; Form2.Image1.Canvas.MoveTo(xn,yn) ; for i := 1 to limiet do begin phi := faktor*(i-1) ; vv := -sin(2*phi) ; xn := xo+Trunc(cos(phi)*sx) ; yn := yo-Trunc(sy*(1+vv)) ; Form2.Image1.Canvas.LineTo(xn,yn) ; end; end; {************************************} procedure MeshPlot; { Plot mesh --------- } var Breed : integer; x, y : integer; m, L : byte; nu, grens : integer; begin Breed := Form2.Image1.Height; Form2.Image1.Canvas.Pen.Color := clBlue; ; grens := (ndr-1)*(nph-1); for nu := 1 to grens do begin PoolWeb(nu, kon[1], kon[2], kon[3], kon[4]); for m := 1 to 4 do begin L := (m mod 4)+1 ; x := Trunc(Breed*0.1*xx(kon[m])); y := Trunc(Breed*0.1*(10-yy(kon[m]))); Form2.Image1.Canvas.MoveTo(x,y); x := Trunc(Breed*0.1*xx(kon[L])); y := Trunc(Breed*0.1*(10-yy(kon[L]))); Form2.Image1.Canvas.LineTo(x,y); end; end; end; {************************************} procedure OutText(schrijf : string); begin with Form2.Image1.Canvas do TextOut(PenPos.x,PenPos.y,schrijf); end; procedure HdBcheck; var { Point on the screen: } waar : vier_ints; xo, yo, sx, sy : integer; {************************************} procedure Cirkel_Numeriek; { Plot calculated velocities on circle ------------------------------------ } var u, v : array of double; ff : double; xn, yn : integer; k : byte; file1 : text; begin SetLength(u,nph-1+1); SetLength(v,nph-1+1); Assign(file1, naam); Reset(file1); for k := 1 to nph-1 do begin Readln(file1, ff); u[k] := ff; Readln(file1, ff); v[k] := ff; end; Close(file1); Form2.Image1.Canvas.Pen.Color := clGreen; { Horizontal component ... } xn := xo+Trunc(sx*0.5*(xx(1)+xx(2))); yn := yo-Trunc(sy*(1+u[1])); Form2.Image1.Canvas.MoveTo(xn,yn) ; for k := 2 to nph-1 do begin xn := xo+Trunc(sx*0.5*(xx(k)+xx(k+1))) ; yn := yo-Trunc(sy*(1+u[k])) ; Form2.Image1.Canvas.LineTo(xn,yn) ; end; { Vertical component ... } xn := xo+Trunc(sx*0.5*(xx(1)+xx(2))); yn := yo-Trunc(sy*(1+v[1])); Form2.Image1.Canvas.MoveTo(xn,yn) ; for k := 2 to nph-1 do begin xn := xo+Trunc(sx*0.5*(xx(k)+xx(k+1))) ; yn := yo-Trunc(sy*(1+v[k])) ; Form2.Image1.Canvas.LineTo(xn,yn) ; end; SetLength(u,0); SetLength(v,0); end; {************************************} begin with Form2.Image1 do begin xo := 50 ; yo := Height-10 ; sx := Width div 4 ; sy := Height div 4 ; end; { Coordinate system: } with Form2.Image1.Canvas do begin Pen.Color := clBlack; MoveTo(xo,yo-sy) ; LineTo(xo+sx,yo-sy) ; MoveTo(xo,yo) ; LineTo(xo,yo-4*sy) ; end; waar.xo := xo ; waar.yo := yo ; waar.sx := sx ; waar.sy := sy ; Writeln('Comparison Exact with Numerical'); Cirkel_Analytisch(waar); Cirkel_Numeriek; with Form2.Image1.Canvas do begin MoveTo(200,48) ; Font.Color := clBlack ; OutText('Comparison between'); Font.Color := clBlue; OutText(' analytical'); Font.Color := clBlack; OutText(' and'); Font.Color := clGreen ; OutText(' numerical'); Font.Color := clBlack ; OutText(' solution'); MoveTo(200,60) ; OutText('for both velocity components at the cylinder surface'); end; end; {************************************} procedure HdBpost1; begin Writeln('The calculated Ideal Flow velocity field'); MeshPlot; VectorPlot(kenmerk+'flow.txt',lengte); with Form2.Image1.Canvas do begin MoveTo(300,30) ; Font.Color := clBlack ; OutText('Velocity field calculated using'); MoveTo(300,42) ; OutText('L.S.FEM with linear '+naam); end; end; procedure HdBpost2; begin Writeln('Potential and Streamlines,'); Writeln('calculated using the Ideal Flow velocity field'); MeshPlot; Form2.Image1.Canvas.Pen.Color := clRed; ; HdBcontour(nivos,kenmerk+'psie.txt'); Form2.Image1.Canvas.Pen.Color := clPurple ; HdBcontour(nivos,kenmerk+'phie.txt'); with Form2.Image1.Canvas do begin MoveTo(300,30) ; Font.Color := clRed ; OutText('Streamlines'); Font.Color := clBlack ; OutText(' and'); Font.Color := clPurple ; OutText(' Potential'); end; end; procedure HdBpost3; begin Writeln('Potential and Streamlines,'); Writeln('calculated using an ab initio procedure'); MeshPlot; Form2.Image1.Canvas.Pen.Color := clPurple ; HdBcontour(nivos,'PhieHdB.txt'); Form2.Image1.Canvas.Pen.Color := clRed ; HdBcontour(nivos,'PsieHdB.txt'); with Form2.Image1.Canvas do begin MoveTo(300,30) ; Font.Color := clRed ; OutText('Streamlines'); Font.Color := clBlack ; OutText(' and'); Font.Color := clPurple ; OutText(' Potential'); end; end; procedure HdBpost4; begin Writeln('Velocity field, derived from'); Writeln('Potential and Streamlines'); MeshPlot; VectorPlot('FlowHdB.txt',lengte); with Form2.Image1.Canvas do begin MoveTo(300,30) ; Font.Color := clBlack ; OutText('Velocity field calculated using'); MoveTo(300,42) ; OutText('both Potential and Streamfunction'); end; end; END.