unit Laplace;

INTERFACE

uses SysUtils;

type
  integers = array of integer;
  vector = array of double;
var
  Nx,Ny : integer;
  b : vector;

procedure Starten;
function nr(i,j : integer) : integer;
procedure Rekenen;
procedure Opschrijven;

IMPLEMENTATION

type
  matrix = array of array of double;
var
  nn,nb1 : integer;
  no : integers;
  S,r : vector;
  E : matrix;

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

procedure four(i,j : integer);
{
  Rectangle of Resistors
}
const
  lambda : double = 1;
var
  dx,dy,Ax,Ay : double;
begin
  SetLength(no,4);
  no[0] := nr(i,j);
  no[1] := nr(i+1,j);
  no[2] := nr(i,j+1);
  no[3] := nr(i+1,j+1);
  dx := 1; dy := 1;
  Ax := lambda/dx*dy/2;
  Ay := lambda/dy*dx/2;
  SetLength(E,4,4);
  E[0,1] := - Ax; E[0,2] := - Ay; E[0,3] := 0; E[0,0] := Ax + Ay;
  E[1,0] := - Ax; E[1,2] := 0; E[1,3] := - Ay; E[1,1] := Ax + Ay;
  E[2,0] := - Ay; E[2,1] := 0; E[2,3] := - Ax; E[2,2] := Ax + Ay;
  E[3,0] := 0; E[3,1] := - Ay; E[3,2] := - Ax; E[3,3] := Ax + Ay;
  SetLength(R,4);
  R[0] := 0; R[1] := 0; R[2] := 0; R[3] := 0;
end;

function lok(i,j : integer) : integer;
begin
  lok := (i-1)*nb1 + j;
end;

procedure Intellen(nds : integer);
{
  Finite Element assembly
  ----------------------- }
var
  m, n : integer;
  ii, jj, ig, jg : integer;
  ij : integer;
begin
{ Use topology: }
  for m := 0 to nds-1 do
  begin
    ig := no[m]+1 ;
  { Half of band: }
    for n := m to nds-1 do
    begin
      jg := no[n]+1 ;
    { Symmetric storage rows: }
      ii := ig ; if jg < ig then ii := jg ;
      if ii > nn then begin
        Writeln('Matrix size exceeded');
      end;
    { Symmetric storage columns: }
      jj := abs(ig-jg)+1 ;
      if jj > nb1 then begin
        Writeln('Bandwidth exceeded');
      end;
    { Account for element: }
      ij := lok(ii,jj);
      S[ij] := S[ij] + E[m,n] ;
    end;
    b[ig] := b[ig] + R[m] ;
  end;
end;

procedure Oplossen;
{
  Solving a linear system of equations
  with positive definite symmetric band matrix
}
var
  i, ii, k, kk, j, jj, n : integer;
  diag, pivot : double;
begin
  if nn < 2 then
    Writeln('Oplossen: nn = '+IntToStr(nn));
  if nb1 < 2 then
    Writeln('Oplossen: nb1 = '+IntToStr(nb1));

  ii := nn-1 ; kk := nb1-1 ;
  for i := 1 to ii do
  begin
    diag := s[lok(i,1)] ;
    if diag = 0 then begin
      Writeln(': 0 on diagonal at: '+IntToStr(i));
    end ;
    if i > nn-nb1 then kk := nn-i ;
    for k := 1 to kk do
    begin
      pivot := s[lok(i,k+1)]/diag ;
      if pivot = 0 then Continue ;
      b[i+k] := b[i+k] - pivot*b[i] ;
      jj := kk-k+1 ;
      for j := 1 to jj do
        s[lok(i+k,j)] := s[lok(i+k,j)] - pivot*s[lok(i,k+j)] ;
    end;
  end;

  kk := nb1-1 ;
  for n := 2 to nn do
  begin
    i := nn-n+2 ;
    b[i] := b[i]/s[lok(i,1)] ;
    if i < nb1 then kk := i-1 ;
    for k := 1 to kk do
      b[i-k] := b[i-k] - b[i]*s[lok(i-k,k+1)] ;
  end;
  b[1] := b[1]/s[lok(1,1)] ;
end;

procedure Schoonmaak;
{
  Clear global matrix & load vector
}
var
  i, j : integer;
begin
  for i := 1 to nn do
  begin
    b[i] := 0;
    for j := 1 to nb1 do
      s[lok(i,j)] := 0;
  end;
end;

procedure Randvoorwaarde;
{
  Insert boundary condition
}
var
  waarde : double;
  ig, i, ib : integer;
begin
  waarde := r[0]/e[0,0] ;
  ig := no[0]+1 ;

  ib := 1 ; if ig > nb1 then ib := ig-nb1+1 ;
  if ig > 1 then
  for i := ib to ig-1 do
  begin
    b[i] := b[i] - s[lok(i,ig-i+1)] * waarde ;
    s[lok(i,ig-i+1)] := 0;
  end;

  ib := nn ; if ig <= nn-nb1 then ib := ig+nb1-1 ;
  if ig < nn then
  for i := ig+1 to ib do
  begin
    b[i] := b[i] - s[lok(ig,i-ig+1)] * waarde ;
    s[lok(ig,i-ig+1)] := 0 ;
  end;

  s[lok(ig,1)] := 1 ;
  b[ig] := waarde ;
end;

procedure Opschrijven;
{
  Write solution to file
}
var
  i,j,k : integer;
  ff : TextFile;
begin
  Assign(ff,'results.txt');
  Rewrite(ff);
  for j := 0 to Ny-1 do
  begin
    for i := 0 to Nx-1 do
    begin
      k := nr(i,j)+1;
      Writeln(ff,i:4,j:4,' ',b[k]) ;
    end;
  end;
  Close(ff);
end;

procedure Rekenen;
{
  Main calculations
}
var
  i,j : integer;
  oo : double;
begin
  Schoonmaak;
{ Bulk elements }
  for j := 0 to Ny-2 do
  begin
    for i := 0 to Nx-2 do
    begin
      four(i,j);
      Intellen(4);
    end;
  end;
{ Boundaries }
  E[0,0] := 1;
  oo := 2*Pi/(Nx-1);
  for i := 0 to Nx-1 do
  begin
    R[0] := cos(oo*i);
    no[0] := nr(i,0);
    Randvoorwaarde;
  end;
  Oplossen;
end;

procedure Starten;
{
  Initialize
}
begin
  nn := Nx*Ny;
  nb1 := Nx+2;
  SetLength(b,nn+1);
  SetLength(s,(nn+1)*(nb1+1));
end;

END.