unit Unit1;
{
  This software has been designed and is CopyLefted
  by Han de Bruijn:
                          (===)
                         @-O^O-@
                          #/_\#
                           ###                    }
interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, Grafisch, Omcirkel, Numeriek, PaulBourke, Contours;

type
  TForm1 = class(TForm)
    Image1: TImage;
    procedure Scheppen(Sender: TObject);
    procedure Toetsdruk(Sender: TObject; var Key: Char);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

var
  FEM : Numerical;
  O : Omtrekken;
  keer : integer;

procedure veld(x,y : double; var u,v : double);
{
  Velocity Field of
  Ideal Flow past a
  circular Cylinder
}
var
  a,n : double;
begin
  a := R0;
  n := sqr(sqr(x)+sqr(y));
  u := 1-sqr(a)*(sqr(x)-sqr(y))/n;
  v := -sqr(a)*2*x*y/n;
end;

procedure Tekenen;
{
  Visualization of
  Finite Element grid
}
const
  s : double = 0.3;
var
  i,j,k,L : integer;
  x,y,u,v : double;
  p,q : punt;

procedure SubRoutine;
begin
  p := q;
  q := plaats(Nummer(i,j));
  x := (p.x+q.x)/2; y := (p.y+q.y)/2;
  k := x2i(x); L := y2j(y);
  veld(x,y,u,v);
  with Form1.Image1.Canvas do
  begin
    Pen.Color := clBlack;
    LineTo(x2i(q.x),y2j(q.y));
    Ellipse(k-2,L-2,k+2,L+2);
    Pen.Color := clRed;
    MoveTo(x2i(x),y2j(y));
    LineTo(x2i(x+s*u),y2j(y+s*v));
    MoveTo(x2i(q.x),y2j(q.y));
  end;
end;

begin
  Form1.Image1.Canvas.Pen.Width := 2;
  for i := 1 to ndr do
  begin
    q := plaats(Nummer(i,1));
    Form1.Image1.Canvas.MoveTo(x2i(q.x),y2j(q.y));
    for j := 2 to nph do
    begin
      SubRoutine;
    end;
  end;
  for j := 1 to nph do
  begin
    q := plaats(Nummer(1,j));
    Form1.Image1.Canvas.MoveTo(x2i(q.x),y2j(q.y));
    for i := 2 to ndr do
    begin
      SubRoutine;
    end;
  end;
{ p := plaats(Nummer(1,nph));
  k := x2i(p.x); L := y2j(p.y);
  with Form1.Image1.Canvas do
  begin
    Pen.Color := clBlack;
    Brush.Color := clBlack;
    Ellipse(k-5,L-5,k+5,L+5);
  end; }
end;

procedure Rekenen;
{
  Calculations
}
{ const
  BIG : double = 1.E+3; }
var
  u,v,x,y : double;
  C1,C2 : punt;
  i,j : integer;
begin
  FEM.nn := 0;
  FEM.nb1 := 0;
  SetLength(FEM.nr,2);
  for j := 1 to nph do
  begin
    for i := 1 to ndr-1 do
    begin
      FEM.nr[0] := Nummer(i,j); 
      FEM.nr[1] := Nummer(i+1,j); 
      FEM.Layout;
    end;
  end;
  for i := 1 to ndr do
  begin
    for j := 1 to nph-1 do
    begin
      FEM.nr[0] := Nummer(i,j); 
      FEM.nr[1] := Nummer(i,j+1); 
      FEM.Layout;
    end;
  end;
  Writeln('Matrix size = ',FEM.nn,' x ',FEM.nb1);
  FEM.Globaal_nul;

  SetLength(FEM.A,2);
  FEM.A[0] := -1; FEM.A[1] := +1;
  for j := 1 to nph do
  begin
    for i := 1 to ndr-1 do
    begin
      C1 := plaats(Nummer(i,j));
      C2 := plaats(Nummer(i+1,j));
      x := (C1.x+C2.x)/2;
      y := (C1.y+C2.y)/2;
      veld(x,y,u,v);
    { \psi_2-\psi_1 = -v(x_2-x_1)+u(y_2-y_1) }
      FEM.rhs := -v*(C2.x-C1.x)+u*(C2.y-C1.y);
      FEM.Element_nul(2);
      FEM.Kleinste_kwadraten(1);
      FEM.nr[0] := Nummer(i,j); 
      FEM.nr[1] := Nummer(i+1,j); 
      FEM.Intellen;
    end;
  end;
  for i := 1 to ndr do
  begin
    for j := 1 to nph-1 do
    begin
      C1 := plaats(Nummer(i,j));
      C2 := plaats(Nummer(i,j+1));
      x := (C1.x+C2.x)/2;
      y := (C1.y+C2.y)/2;
      veld(x,y,u,v);
    { \psi_2-\psi_1 = -v(x_2-x_1)+u(y_2-y_1) }
      FEM.rhs := -v*(C2.x-C1.x)+u*(C2.y-C1.y);
      FEM.Element_nul(2);
      FEM.Kleinste_kwadraten(1);
      FEM.nr[0] := Nummer(i,j); 
      FEM.nr[1] := Nummer(i,j+1); 
      FEM.Intellen;
    end;
  end;
{ i := Nummer(1,nph);
  k := FEM.lok(i,i);
  FEM.s[k-1] := BIG; }
  FEM.nr[0] := Nummer(1,nph);
  FEM.e[0,0] := 1; FEM.r[0] := 0;
  FEM.Randvoorwaarde;
  FEM.Oplossen;
end;

procedure TForm1.Scheppen(Sender: TObject);
{
  On Create
}
begin
  TV(Form1.Image1);
  ClearDevice;
  xmin := -0.3; xmax := 5.3;
  ymin := -0.3; ymax := 5.3;
  ndr := 10; nph := 10;
  R0 := 1; R1 := 5;

  Rekenen;

  O := Omtrekken.Create;
  O.afmetingen(Wijd,Hoog);
  O.Scherm := Form1.Image1.Picture.Bitmap;
  keer := 0;
end;

procedure Rasteren(q : isoquad);
{
  Show Grid with Contours
}
var
  xm,ym : double;
begin
  xm := (q[1].x+q[2].x+q[3].x+q[4].x)/4;
  ym := (q[1].y+q[2].y+q[3].y+q[4].y)/4;
  with Form1.Image1.Canvas do
  begin
    Pen.Width := 1;
    Pen.Color := clGreen;
    MoveTo(x2i(q[1].x),y2j(q[1].y));
    LineTo(x2i(q[2].x),y2j(q[2].y));
    LineTo(x2i(q[4].x),y2j(q[4].y));
    LineTo(x2i(q[3].x),y2j(q[3].y));
    LineTo(x2i(q[1].x),y2j(q[1].y));
    MoveTo(x2i(q[1].x),y2j(q[1].y));
    LineTo(x2i(xm),y2j(ym));
    MoveTo(x2i(q[2].x),y2j(q[2].y));
    LineTo(x2i(xm),y2j(ym));
    MoveTo(x2i(q[3].x),y2j(q[3].y));
    LineTo(x2i(xm),y2j(ym));
    MoveTo(x2i(q[4].x),y2j(q[4].y));
    LineTo(x2i(xm),y2j(ym));
    Pen.Width := 2;
    Pen.Color := clBlack;
  end;
end;

function psi(x,y : double) : double;
begin
  psi := (1-sqr(R0)/(sqr(x)+sqr(y)))*y;
end;

procedure Isolines;
{
  Producing IsoLines of
  given Velocity Field
}
const
  meest : integer = 40;
  BIG : double = 1E6;
var
  i,j,k,n : integer;
  C : array[1..4] of punt;
  q : isoquad;
  bmin,bmax : double;
  p : punt;
  x,y,nivo,max : double;
begin
  bmin := BIG; bmax := 0;
  for k := 0 to FEM.nn-1 do
  begin
  { Writeln(FEM.b[k]); }
    if FEM.b[k] < bmin then bmin := FEM.b[k];
    if FEM.b[k] > bmax then bmax := FEM.b[k];
  end;
{ Writeln(bmin,bmax); }
  for k := 0 to FEM.nn-1 do
  begin
    FEM.b[k] := (FEM.b[k]-bmin)/(bmax-bmin);
  end;
{ Writeln(bmin,bmax); }
  Setup_Contours;
  Beeld := Form1.Image1;
  for i := 1 to ndr-1 do
  begin
    for j := 1 to nph-1 do
    begin
      C[1] := plaats(Nummer(i,j));
      C[2] := plaats(Nummer(i+1,j));
      C[3] := plaats(Nummer(i,j+1));
      C[4] := plaats(Nummer(i+1,j+1));
      for k := 1 to 4 do
      begin
        q[k].x := C[k].x;
        q[k].y := C[k].y;
      end;
      q[1].z := FEM.b[Nummer(i,j)-1];
      q[2].z := FEM.b[Nummer(i+1,j)-1];
      q[3].z := FEM.b[Nummer(i,j+1)-1];
      q[4].z := FEM.b[Nummer(i+1,j+1)-1];
      Rasteren(q);
      Contour_Vierhoek(meest,q);
    end;
  end;
{ Zero Streamline }
  p := plaats(Nummer(1,nph));
  Form1.Image1.Canvas.MoveTo(x2i(p.x),y2j(p.y));
  for j := nph-1 downto 1 do
  begin
    p := plaats(Nummer(1,j));
    Form1.Image1.Canvas.LineTo(x2i(p.x),y2j(p.y));
  end;
  for i := 2 to ndr do
  begin
    p := plaats(Nummer(i,1));
    Form1.Image1.Canvas.LineTo(x2i(p.x),y2j(p.y));
  end;
  if keer = 1 then Exit;
{ Exact solution contours }
  max := psi(0,R1);
  for n := 0 to meest do
  begin
    nivo := n/meest*max;
    for j := 0 to Hoog-1 do
    begin
      for i := 0 to Wijd-1 do
      begin
        x := i2x(i); y := j2y(j);
        O.funktie[i,j] := 0;
	if sqr(x)+sqr(y) < sqr(R0) then Continue;
	if sqr(x)+sqr(y) > sqr(R1) then Continue;
	if x < 0 then Continue;
        O.funktie[i,j] := psi(x,y)-nivo;
      end;
    end;
    O.Rondom(false);
    for k := 0 to O.Aantal do
      O.Schetsen(k,clRed,2);
  end;
end;

procedure TForm1.Toetsdruk(Sender: TObject; var Key: Char);
{
  On KeyPress
}
begin
  ClearDevice;
  Case keer of
  0 : Tekenen;
  1 : Isolines;
  2 : Isolines;
  end;
  keer := keer + 1;
  if keer = 3 then keer := 0;
end;

end.