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

Uses Numeriek, Generaal;

const
  NDX = 18 ; NDY = 3;
var
  FEM : Numerical;
  { numerical method }
  crd : array of array of punt;
  { Quadrilateral coordinates } 

procedure Rekenen; { Calculations }
procedure topology(i,j : integer);

IMPLEMENTATION

function plaats(i,j : integer) : punt;
{
  Global (x,y) Coordinates
}
begin
  plaats := crd[i-1,j-1]; 
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);
{ Writeln((x2-x1)*(y4-y3)-(x4-x3)*(y2-y1)); }
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_boven(i : integer);
{ 
  Impermeable: x.u + y.v = 0
}
var
  no : integer;
  p,q : punt;
  x21,y21 : double;
begin
  FEM.Element_nul(2);
  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;
  FEM.a[0] := y21 ; FEM.a[1] := -x21;
  FEM.rhs := 0;
  FEM.Kleinste_kwadraten(1);
  FEM.Intellen;

  FEM.Element_nul(2);
  no := (2*i-1)*NDY+i+NDY;
  FEM.nr[0] := no*2-1;
  FEM.nr[1] := no*2;
  p := plaats(i,NDY+1);
  q := plaats(i+1,NDY+1);
  x21 := q.x-p.x; y21 := q.y-p.y;
  FEM.a[0] := y21 ; FEM.a[1] := -x21;
  FEM.rhs := 0;
  FEM.Kleinste_kwadraten(1);
  FEM.Intellen;
end;

procedure links(j : integer);
{ 
  Inflow: (u,v)_n = given
}
var
  no : integer;
  p,q : punt;
  x21,y21 : double;
begin
  FEM.Element_nul(2);
  no := j;
  FEM.nr[0] := no*2-1;
  FEM.nr[1] := no*2;
  p := plaats(1,j);
  q := plaats(1,j+1);
  x21 := q.x-p.x; y21 := q.y-p.y;
  FEM.a[0] := y21 ; FEM.a[1] := -x21;
  FEM.rhs := -sqrt(sqr(x21)+sqr(y21));
  FEM.Kleinste_kwadraten(1);
  FEM.Intellen;
end;

procedure rechts(j : integer);
{ 
  Outflow: (u,v) // normal
}
var
  no : integer;
  p,q : punt;
  x21,y21 : double;
begin
  FEM.Element_nul(2);
  no := NDX*(2*NDY+1)+j;
  FEM.nr[0] := no*2-1;
  FEM.nr[1] := no*2;
  p := plaats(NDX+1,j);
  q := plaats(NDX+1,j+1);
  x21 := q.x-p.x; y21 := q.y-p.y;
  FEM.a[0] := x21 ; FEM.a[1] := y21;
  FEM.rhs := 0;
  FEM.Kleinste_kwadraten(1);
  FEM.Intellen;
end;

procedure Rekenen;
{
  Calculations
}
var
  i,j,k : integer;
  uit : TextFile;
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);
  SetLength(FEM.a,2);
  for i := 1 to NDX do
    onder_boven(i);
  for j := 1 to NDY do
    links(j);
  for j := 1 to NDY do
    rechts(j);

  FEM.Oplossen;

  AssignFile(uit,'DuctFlow.txt');
  Rewrite(uit);
  for k := 0 to FEM.nn-1 do
  begin
    Writeln(uit,FEM.b[k]);
  end;
  CloseFile(uit);
end;

END.