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

Uses Numeriek;

type
  punt = record
{ point in space }
  x,y : double;
end;

var
  NDX,NDY : integer;
  { 2-D grid size }
  FEM : Numerical;
  { numerical method }
  vgl : array of double;
  { analytical solution }

procedure Rekenen; { Calculations }

IMPLEMENTATION

const
  R0 : double = 1;
  { inner radius }
  R1 : double = 10;
  { outer radius }

function plaats(i,j : integer) : punt;
{
  Global (x,y) Coordinates
}
var
  R,phi : double;
  p : punt;
begin
  R := R0*exp((j-1)/NDY*ln(R1/R0));
  phi := 0.5*Pi*(i-1)/NDX;
  p.x := R*cos(phi) ;
  p.y := R*sin(phi) ;
  plaats := p;
end;

procedure element(i,j : integer);
{
  Bulk Element for Ideal Flow
}
var
  A,B,C,D : punt;
  x1,x2,x3,x4 : double;
  y1,y2,y3,y4 : double;
begin
  A := plaats(i,j); B := plaats(i+1,j);
  C := plaats(i,j+1); D := plaats(i+1,j+1);

  x1 := (A.x+C.x)/2; y1 := (A.y+C.y)/2;
  x2 := (B.x+D.x)/2; y2 := (B.y+D.y)/2;
  x3 := (A.x+B.x)/2; y3 := (A.y+B.y)/2;
  x4 := (C.x+D.x)/2; y4 := (C.y+D.y)/2;

  FEM.Element_nul(8);
  FEM.rhs := 0;
  FEM.a[0] := -(y4-y3) ; FEM.a[1] := +(x4-x3);
  FEM.a[2] := +(y4-y3) ; FEM.a[3] := -(x4-x3);
  FEM.a[4] := +(y2-y1) ; FEM.a[5] := -(x2-x1);
  FEM.a[6] := -(y2-y1) ; FEM.a[7] := +(x2-x1);
  FEM.Kleinste_kwadraten(1);
  FEM.a[0] := -(x4-x3) ; FEM.a[1] := -(y4-y3);
  FEM.a[2] := +(x4-x3) ; FEM.a[3] := +(y4-y3);
  FEM.a[4] := +(x2-x1) ; FEM.a[5] := +(y2-y1);
  FEM.a[6] := -(x2-x1) ; FEM.a[7] := -(y2-y1);
  FEM.Kleinste_kwadraten(1);
  FEM.a[0] := +1; FEM.a[1] := 0; FEM.a[2] := +1; FEM.a[3] := 0;
  FEM.a[4] := -1; FEM.a[5] := 0; FEM.a[6] := -1; FEM.a[7] := 0;
  FEM.Kleinste_kwadraten(1);
  FEM.a[0] := 0; FEM.a[1] := +1; FEM.a[2] := 0; FEM.a[3] := +1;
  FEM.a[4] := 0; FEM.a[5] := -1; FEM.a[6] := 0; FEM.a[7] := -1;
  FEM.Kleinste_kwadraten(1);
end;

procedure topology(i,j : integer);
var
  k : integer;
  no : array[1..4] of integer;
begin
  no[1] := (i-1)*(2*NDY+1)+j;
  no[2] := i*(2*NDY+1)+j;
  no[3] := (2*i-1)*NDY+j+(i-1);
  no[4] := (2*i-1)*NDY+j+i;
{ for k := 1 to 4 do
    Write(no[k]:4);
  Writeln; }
  for k := 1 to 4 do
  begin
    FEM.nr[2*(k-1)] := no[k]*2-1;
    FEM.nr[2*(k-1)+1] := no[k]*2;
  end;
end;

procedure onder(i : integer);
{ 
  Circle: x.u + y.v = 0
}
var
  no : integer;
  p,q : punt;
  x21,y21,noemer : double;
begin
  no := (2*i-1)*NDY+i;
  FEM.nr[0] := no*2-1;
  FEM.nr[1] := no*2;
  p := plaats(i,1);
  q := plaats(i+1,1);
  x21 := q.x-p.x; y21 := q.y-p.y;
  noemer := 1/(x21*x21 + y21*y21);
  FEM.e[0,0] := y21*y21*noemer;
  FEM.e[1,1] := x21*x21*noemer ;
  FEM.e[0,1] := -x21*y21*noemer;
  FEM.e[1,0] := -x21*y21*noemer ;
  FEM.Intellen;
end;

procedure boven(i : integer);
{ 
  Outflow: u = 1
}
var
  no : integer;
begin
  no := (2*i-1)*NDY+i+NDY;
  FEM.nr[0] := no*2-1;
  FEM.nr[1] := no*2;
  FEM.e[0,0] := 1;
  FEM.r[0] := 1;
  FEM.Intellen;
end;

procedure links(j : integer);
{ 
  Inflow: v = 0
}
var
  no : integer;
begin
  no := j;
  FEM.nr[0] := no*2-1;
  FEM.nr[1] := no*2;
  FEM.e[1,1] := 1;
  FEM.Intellen;
end;

procedure rechts(j : integer);
{ 
  Symmetry: v = 0
}
var
  no : integer;
begin
  no := NDX*(2*NDY+1)+j;
  FEM.nr[0] := no*2-1;
  FEM.nr[1] := no*2;
  FEM.e[1,1] := 1;
  FEM.Intellen;
end;

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

function half(C1,C2 : punt) : punt;
var
  C : punt;
begin
  C.x := (C1.x+C2.x)/2;
  C.y := (C1.y+C2.y)/2;
  half := C;
end;

procedure Vergelijken;
{
  Analytical solution
  at same nodal points
}
var
  i,j,k : integer;
  A,B,C,D : punt;
  p : array[1..4] of punt;
  w : punt;
begin
  SetLength(vgl,FEM.nn);
  SetLength(FEM.nr,8);
  for i := 1 to NDX do
  begin
    for j := 1 to NDY do
    begin
      A := plaats(i,j); B := plaats(i+1,j);
      C := plaats(i,j+1); D := plaats(i+1,j+1);
      p[1] := half(A,C); p[2] := half(B,D);
      p[3] := half(A,B); p[4] := half(C,D);
      topology(i,j);
      for k := 1 to 4 do
      begin
        veld(p[k],w);
	vgl[FEM.nr[2*(k-1)]-1] := w.x;
	vgl[FEM.nr[2*(k-1)+1]-1] := w.y;
      end;
    end;
  end;
end;

procedure Rekenen;
{
  Calculations
}
var
  i,j,k : integer;
begin
  FEM.nn := 0;
  FEM.nb1 := 0;
  SetLength(FEM.nr,8);
  for i := 1 to NDX do
  begin
    for j := 1 to NDY do
    begin
      topology(i,j);
      FEM.Layout;
    end;
  end;
  Writeln('Matrix size = ',FEM.nn,' x ',FEM.nb1);
  FEM.Globaal_nul;
 
  SetLength(FEM.A,8);
  for i := 1 to NDX do
  begin
    for j := 1 to NDY do
    begin
      element(i,j);
      topology(i,j);
      FEM.Intellen;
    end;
  end;

  SetLength(FEM.nr,2);
  FEM.Element_nul(2);
  for i := 1 to NDX do
    onder(i);
  FEM.Element_nul(2);
  for i := 1 to NDX do
    boven(i);
  FEM.Element_nul(2);
  for j := 1 to NDY do
    links(j);
  FEM.Element_nul(2);
  for j := 1 to NDY do
    rechts(j);

  FEM.Oplossen;
  Vergelijken;

  Writeln('   Numerical            =    Analytical');
  for k := 0 to FEM.nn-1 do
  begin
    Writeln(FEM.b[k],' = ',vgl[k]);
  end;
end;

END.