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

  All I ask is to be credited when it is appropriate.
}
interface

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

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

var
  Form1: TForm1;

implementation

{$R *.dfm}

type
  vector = array of double;
  matrix = array of array of double;
const
  grid : integer = 10;
var
  O : Omtrekken;
  b : vector;
  max : double;

procedure Oplossen(A : matrix; var b : vector);
{
  Simple Gauss Solver
}
var
  N,i,j,k : integer;
  p,s : double;
begin
  N := Length(b);
{ Elimination }
  for k := 0 to N-1 do
  begin
    p := A[k,k];
    for i := k+1 to N-1 do
    begin
      if A[i,k] = 0 then Continue;
      A[i,k] := A[i,k]/p;
      for j := k+1 to N-1 do
        A[i,j] := A[i,j] - A[i,k]*A[k,j];
      b[i] := b[i] - A[i,k]*b[k];
    end;
  end;
{ Backsubstitution }
  for k := N-1 downto 0 do
  begin
    s := b[k];
    for j := k+1 to N-1 do
      s := s - A[k,j]*b[j];
    b[k] := s/A[k,k];
  end;
end;

function nr(i,j : integer) : integer;
begin
  nr := grid*j+i;
end;

procedure Berekenen;
{
   Calculations
}
const
  BIG : double = 1.E+6;
var
  k,i,j,x,y : integer;
  M,N : integer;
  u,v : double;
  A,S : matrix;
  r : vector;
begin
  N := grid*grid;
  M := 2*(grid-1)*grid;
  SetLength(A,M,N);
  SetLength(r,M);
{ Overdetermined system }
  for k := 0 to M-1 do
  for i := 0 to N-1 do
    A[k,i] := 0;
  for k := 0 to M-1 do
    r[k] := 0;
  k := 0;
  for y := 0 to grid-1 do
  begin
    for x := 0 to grid-2 do
    begin
    { d\psi/dx = -v }
      v := -y;
      A[k,nr(x,y)] := -1;
      A[k,nr(x+1,y)] := +1;
      r[k] := -v; 
      k := k+1;
    end;
  end;
  for x := 0 to grid-1 do
  begin
    for y := 0 to grid-2 do
    begin
    { d\psi/dy = u }
      u := x;
      A[k,nr(x,y)] := -1;
      A[k,nr(x,y+1)] := +1;
      r[k] := u; 
      k := k+1;
    end;
  end;
{ Least Squares FEM }
  SetLength(S,N,N);
  SetLength(b,N);
  for i := 0 to N-1 do
  for j := 0 to N-1 do
    S[i,j] := 0;
  for i := 0 to N-1 do
  for j := 0 to N-1 do
  for k := 0 to M-1 do
    S[i,j] := S[i,j] + A[k,i]*A[k,j];
  for i := 0 to N-1 do
    b[i] := 0;
  for i := 0 to N-1 do
  for k := 0 to M-1 do
    b[i] := b[i] + A[k,i]*r[k];
  i := grid*(grid-1);
  S[i,i] := BIG;
  Oplossen(S,b);
  max := 0;
  for i := 0 to N-1 do
    if b[i] > max then max := b[i];
end;

procedure Tekenen;
{
  Contour plot
}
const
  veel : integer = 25;
var
  i,j,L,U,n,k : integer;
  x,y,f0,f1,f2,f3,dx,dy,xi,eta,f,nivo : double;
begin
  dx := 1; dy := 1;
  for n := 0 to veel do
  begin
    nivo := n/veel*max;
    for j := 0 to Hoog-1 do
    begin
      for i := 0 to Wijd-1 do
      begin
        x := i2x(i); y := j2y(j);
        L := Trunc(x); U := Trunc(y);
        if L > grid-2 then L := grid-2;
        if U > grid-2 then U := grid-2;
      { Bilinear interpolation }
        f0 := b[nr(L,U)];
        f1 := b[nr(L+1,U)];
        f2 := b[nr(L,U+1)];
        f3 := b[nr(L+1,U+1)];
        xi := (x-L)/dx; eta := (y-U)/dy;
        f := (1-xi)*(1-eta)*f0 + xi*(1-eta)*f1
           + (1-xi)*eta*f2 + xi*eta*f3;
        O.funktie[i,j] := f-nivo;
      end;
    end;
    O.Rondom(false);
    for k := 0 to O.Aantal do
      O.Schetsen(k,clBlack,2);
  end;
end;

procedure TForm1.Toetsdruk(Sender: TObject; var Key: Char);
{
  On KeyPress
}
begin
  Berekenen;
  ClearDevice;
  Tekenen;
end;

procedure Raster;
{
  Grid plot
}
var
  k : integer;
begin
  for k := 0 to grid-1 do
  with Form1.Image1.Canvas do
  begin
    Pen.Color := clRed;
    MoveTo(x2i(0),y2j(k));
    LineTo(x2i(grid-1),y2j(k));
  end;
  for k := 0 to grid-1 do
  with Form1.Image1.Canvas do
  begin
    Pen.Color := clGreen;
    MoveTo(x2i(k),y2j(0));
    LineTo(x2i(k),y2j(grid-1));
  end;
end;

procedure TForm1.Scheppen(Sender: TObject);
{
  On Create
}
begin
  xmin := -0.1; xmax := grid-1+0.1;
  ymin := -0.1; ymax := grid-1+0.1;
  TV(Form1.Image1);
  O := Omtrekken.Create;
  ClearDevice;
  O.afmetingen(Wijd,Hoog);
  O.Scherm := Form1.Image1.Picture.Bitmap;
  Form1.Image1.Canvas.Pen.Width := 2;
  Raster;
end;

end.