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

Uses SysUtils;

type
  Numerical = object
  public
    a : array of double;
    rhs : double;
    e : array of array of double;
    r : array of double;
    nr : array of integer; { = 1,2,3, .. }
    b : array of double;
    nn : integer ; nb1 : integer ;
    s : array of double;

    procedure Element_nul(nds : integer);
    procedure Print_Element;
    procedure Kleinste_Kwadraten(gewicht : double);
    function Restant : double;

    procedure Layout;
    procedure Globaal_nul;
    procedure Intellen;
    procedure Oplossen;
    procedure Randvoorwaarde;
    function lok(i,j : integer) : integer;
end;

IMPLEMENTATION

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

procedure Numerical.Print_Element;
var
  i,j,nds : integer;
begin
  nds := Length(e);
  for i := 0 to nds-1 do
  begin
    for j := 0 to nds-1 do
      Write(e[i][j]:7:2,' ');
    Writeln(' ',r[i]:7:2);
  end;
  Writeln;
end;

procedure Numerical.Kleinste_Kwadraten;
{
  Least Squares Finite Element Method
  ----------------------------------- }
var
  sum : double;
  i,j,nds : integer;
begin
  nds := Length(a);
{ Trace weighting: }
  sum := 0 ;
  for i := 0 to nds-1 do
    sum := sum + a[i]*a[i] ;
  sum := gewicht*gewicht/sum;

{ Add to element: }
  for i := 0 to nds-1 do
  begin
    r[i] := r[i] + a[i]*rhs*sum;
    for j := 0 to nds-1 do
      e[i][j] := e[i][j] + a[i]*a[j]*sum;
  end;
end;

procedure Numerical.Globaal_nul;
var
  i, j : integer;
begin
  SetLength(b,nn);
  SetLength(s,nn*nb1);
{ Clear global matrix and vector: }
  for i := 1 to nn do
  begin
    b[i-1] := 0;
    for j := 1 to nb1 do
      s[lok(i,j)-1] := 0;
  end;
end;

procedure Numerical.Element_nul;
var
  m, n : integer;
begin
  SetLength(r,nds);
  SetLength(e,nds,nds);
{ Clear element: }
  for m := 0 to nds-1 do
  begin
    r[m] := 0;
    for n := 0 to nds-1 do
      e[m][n] := 0 ;
  end;
end;

function Numerical.Restant;
var
  i,j,nds : integer;
  rest : double;
begin
  rest := 0;
  nds := Length(e);
  for i := 0 to nds-1 do
  for j := 0 to nds-1 do
  begin
    rest := rest + b[nr[i]-1] * e[i][j] * b[nr[j]-1];
  end;
  Restant := rest;
end;

procedure Numerical.Layout;
{
  System Matrix Layout
  -------------------- }
var
  m,n,nds : integer;
  ii,jj,ig,jg : integer;
begin
{ Use topology: }
  nds := Length(nr);
  for m := 0 to nds-1 do
  begin
    ig := nr[m] ;
  { Half of band: }
    for n := m to nds-1 do
    begin
      jg := nr[n] ;
    { Symmetric storage rows: }
      ii := ig ; if jg < ig then ii := jg ;
      if ii > nn then nn := ii;
    { Symmetric storage columns: }
      jj := abs(ig-jg)+1 ;
      if jj > nb1 then nb1 := jj;
    end;
  end;
end;

procedure Numerical.Intellen;
{
  Finite Element assembly
  ----------------------- }
var
  m,n,nds : integer;
  ii,jj,ig,jg : integer;
  ij : integer;
begin
{ Use topology: }
  nds := Length(e);
  for m := 0 to nds-1 do
  begin
    ig := nr[m] ;
  { Half of band: }
    for n := m to nds-1 do
    begin
      jg := nr[n] ;
      ii := ig ; if jg < ig then ii := jg ;
      jj := abs(ig-jg)+1 ;
    { Account for element: }
      ij := lok(ii,jj)-1;
      s[ij] := s[ij] + e[m][n] ;
    end;
    b[ig-1] := b[ig-1] + r[m] ;
  end;
end;

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

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

procedure Numerical.Randvoorwaarde;
var
  waarde : double;
  ig,i,ib : integer;
begin
  waarde := r[0]/e[0,0] ; ig := nr[0] ;
  ib := 1 ; if ig > nb1 then ib := ig-nb1+1 ;
  if ig > 1 then
  for i := ib to ig-1 do
  begin
    b[i-1] := b[i-1] - s[lok(i,ig-i+1)-1] * waarde ;
    s[lok(i,ig-i+1)-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-1] := b[i-1] - s[lok(ig,i-ig+1)-1] * waarde ;
    s[lok(ig,i-ig+1)-1] := 0 ;
  end;
  s[lok(ig,1)-1] := 1 ;
  b[ig-1] := waarde ;
end;

END.