Unit OmCirkel;

INTERFACE

type
  punt = record
  x,y : double;
  end;
var
{ # GridLines in Polar Coordinates }
  ndr,nph : integer; 
{ Inner / Outer Mesh Radius }
  R0,R1 : double;

function plaats(k : integer) : punt;
function Nummer(i,j : integer) : integer;
procedure Ontbind(nu : integer ; var i,j : integer);
procedure Knopen_Midden(nr1,nr2 : integer ; var midden : integer) ;
procedure Midden_Knopen(midden : integer ; var nr1,nr2 : integer) ;
procedure PoolWeb(nu : integer ; var n1,n2,n3,n4 : integer);

IMPLEMENTATION

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

function Nummer(i,j : integer) : integer;
{
  F.D. grid numbers -> F.E. global number
}
begin
  Nummer := (i-1)*nph+j ;
end;

procedure Ontbind(nu : integer ; var i,j : integer);
{
  F.E. global number -> F.D. grid numbers
}
begin
  i := ((nu-1) div nph)+1 ;
  j := ((nu-1) mod nph)+1 ;
end;

procedure Knopen_Midden(nr1,nr2 : integer ; var midden : integer) ;
{
  Node numbers -> Midside number
}
var
  i1,j1,i2,j2 : integer;
  goed : boolean;

procedure Orde(var p,q : integer);
var
  h : integer;
begin
  if p < q then Exit;
  h := p ; p := q ; q := h ;
end;

begin
  midden := 0 ;
  Orde(nr1,nr2) ;
  Ontbind(nr1, i1,j1) ; Ontbind(nr2, i2,j2) ;
  goed := (i1 = i2) xor (j1 = j2) ;
  if (not goed) then Exit ;
  if i1 = i2 then midden := (2*nph-1)*(i1-1)+j1 ;
  if j1 = j2 then midden := (2*nph-1)*i1+j1-nph ;
end;

procedure Midden_Knopen(midden : integer ; var nr1,nr2 : integer) ;
{
  Midside number -> Node numbers
}
var
  x,y,i1,j1,i2,j2 : integer;

begin
  x := ((midden-1) div (2*nph-1)) + 1 ;
  y := ((midden-1) mod (2*nph-1)) + 1 ;
  if y < nph then
  begin
    j1 := y ; j2 := y + 1 ;
    i1 := x ; i2 := x ;
  end else begin
    j1 := y - nph + 1 ; j2 := j1 ;
    i1 := x ; i2 := i1 + 1 ;
  end;
  nr1 := Nummer(i1,j1);
  nr2 := Nummer(i2,j2);
end;

procedure PoolWeb(nu : integer ; var n1,n2,n3,n4 : integer);
{
  Polar mesh topology
  ------------------- }
var
  i, j : integer;
  OK : boolean;
begin
  i := ((nu-1) div (nph-1))+1 ;
  j := ((nu-1) mod (nph-1))+1 ;
  n1 := 0 ; n2 := 0 ; n3 := 0 ; n4 := 0 ;
  OK := (i > 0) and (i < ndr) and (j > 0) and (j < nph);
  if (not OK) then Exit;
  n2 := (i-1)*nph+j ; n1 :=(i-1)*nph+j+1 ;
  n3 := i*nph+j ; n4 :=i*nph+j+1 ;
end;

END.