-- pde2sin_la.adb  Galerkin FEM
-- solve ux(x,y) + uy(x,y) = c(x,y)
-- boundary conditions computed using u(x)
-- analytic solution may be given by u(x) = sin(x-y)
-- Gauss-Legendre integration used
-- solve for Uij numerical approximation U(x_i,y_j) of u(x_i,y_j)
-- k * U = f, solve simultaneous equations for U

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with Calendar;
use  Calendar;
with Real_Arrays; -- types real, real_vector, real_matrix
use  Real_Arrays;
with laphi;       -- phi, phip, phipp, phippp, phipppp
use  laphi;
with simeq;       -- function simeq
with gaulegf;     -- function gaulegf

procedure pde2sin_la is
  nx : Integer := 19;
  ny : Integer := 18;
  kg : Real_Matrix(1..nx*ny,1..nx*ny); --((i-1)*ny+ii,(j-1)*ny+jj)
  xg : Real_Vector(1..nx);
  yg : Real_Vector(1..ny);
  fg : Real_Vector(1..nx*ny); -- (i-1)*ny+ii
  ug : Real_Vector(1..nx*ny);
  Ua : Real_Vector(1..nx*ny);
  Pi : Real := 3.1415926535897932384626433832795028841971;
  xmin : real := 0.0;
  xmax : real := 4.0*Pi;
  ymin : real := 0.0;
  ymax : real := 4.0*Pi;
  x, y, hx, hy : real;
  npx , npy : Integer := 23;
  xx : Real_Vector(1..npx); -- for Gauss-Legendre
  wx : Real_Vector(1..npx);
  yy : Real_Vector(1..npy); -- for Gauss-Legendre
  wy : Real_Vector(1..npy);
  i, ii : Integer;
  val, err, avgerr, maxerr : real;
  time_start, time_now : duration;

  function FF(x,y:real) return real is -- forcing function
  begin
    return 0.0;
  end FF;

  function uana(x,y:real) return real is -- analytic solution and boundary
  begin
    return Sin(X-Y);
  end uana;

  function galk(x,y:Real; i,j,ii,jj:integer) return real is
                                       -- Galerkin k stiffness function
  begin                                -- for this specific problem
    return (  phip(x,j,nx,xg)*  phi(y,jj,ny,yg) +
              phi(x,j,nx,xg) * phip(y,jj,ny,yg)  ) *
           phi(x,i,nx,xg)*  phi(y,ii,ny,yg);
  end galk;

  function galf(x,y:Real; i,ii:integer) return real is -- Galerkin f function
  begin
    return FF(x,y)*phi(x,i,nx,xg)*phi(y,ii,ny,yg);
  end galf;

  package Real_IO is new Float_IO(Real);
  use Real_IO;

  -- write_soln2  for gnuplot
  procedure Write_Soln2(U : Real_Vector; File_Name : String) is
    Fout : File_Type;
  begin
    Create(Fout, Out_File, File_Name);

    Put_Line("writing "&File_name);
    for  i in 1..nx loop
      for ii in 1..ny loop
        Put(Fout," ");
        Put(Fout,xg(i),2,6,0);
        Put(Fout," ");
        Put(Fout,yg(ii),2,6,0);
        Put(Fout," ");
        Put(Fout,U((i-1)*ny+ii),4,6,0);
        New_Line(Fout);
      end loop; -- ii
      New_Line(Fout);
    end loop; -- i
    Close(Fout);
    Put_Line("finished writing "&File_Name);
  end Write_Soln2;

begin
  Put_Line("pde2sin_la.adb running ");
  Put_Line("solve  ux(x,y)  +  uy(x,y)  = c(x,y)");
  Put_Line("c(x,y) = 0.0");
  Put_Line("boundary conditions computed using u(x)");
  Put_Line("analytic solution may be given by u(x) = sin(x-y)");
  Put_Line("Gauss-Legendre integration used");
  New_Line;
  Put("xmin="); Put(xmin,2,2,0);
  Put(", xmax="); Put(xmax,2,2,0);
  Put(", ymin="); Put(ymin,2,2,0);
  Put(", ymax="); Put(ymax,2,2,0);
  New_Line;
  Put_Line("nx="&Integer'Image(nx)&", ny="&Integer'Image(ny));
  time_start := calendar.seconds(calendar.clock);

  -- initialize xg array. does not have to be uniform spacing
  hx := (xmax-xmin)/real(nx-1);
  Put_Line("x grid and analytic solution at ymin ");
  for i in 1..nx loop
    xg(i) := xmin+Real(i-1)*hx;
    Ua(i) := uana(xg(i),ymin);
    Put("i="&Integer'Image(i)&", Ua(");
    Put(xg(i),2,5,0);
    Put(",");
    Put(ymin,2,5,0);
    Put(")=");
    Put(Ua(i),3,5,0);
    New_Line;
  end loop;
  New_Line;

  -- initialize yg array. does not have to be uniform spacing
  hy := (ymax-ymin)/real(ny-1);
  Put_Line("y grid and analytic solution at xmin");
  for ii in 1..ny loop
    yg(ii) := ymin+Real(ii-1)*hy;
    Ua(ii) := uana(xmin,yg(ii));
    Put("ii="&Integer'Image(ii)&", Ua(");
    Put(xmin,2,5,0);
    Put(",");
    Put(yg(ii),2,5,0);
    Put(")=");
    Put(Ua(ii),3,5,0);
    New_Line;
  end loop;
  New_Line;

  -- put solution in Ua vector
  for i in 1..nx loop
    x := xmin+Real(i-1)*hx;
    for ii in 1..ny loop
      y := ymin+Real(ii-1)*hy;
      Ua((i-1)*ny+ii) := uana(x,y);
      -- Put("solution at i="&Integer'Image(I));
      -- Put(", x="); Put(X,2,3,0);
      -- Put(", ii="&Integer'Image(Ii));
      -- Put(", y="); Put(Y,2,3,0);
      -- Put(" is "); Put(Ua((i-1)*ny+ii),3,3,0);
      -- New_Line;
    end loop;
  end loop;

  -- initialize k and f
  for i in 1..nx loop
    for j in 1..nx loop
      for ii in 1..ny loop
        for jj in 1..ny loop
          kg((i-1)*ny+ii,(j-1)*ny+jj) := 0.0;
        end loop;
      end loop;
    end loop;
  end loop;
  for i in 1..nx loop
    for ii in 1..ny loop
      fg((i-1)*ny+ii) := 0.0;
    end loop;
  end loop;

  -- set boundary conditions, four sides
  for i in 1..nx loop
    x := xg(i);
    ii := 1;
    y := yg(ii);
    kg((i-1)*ny+ii,(i-1)*ny+ii) := 1.0;
    fg((i-1)*ny+ii) := uana(x,y);
    -- Put("boundary i="&Integer'Image(i));
    -- Put(", x="); Put(x,2,3,0);
    -- Put(", ii="&Integer'Image(ii));
    -- Put(", y="); Put(y,2,3,0);
    -- Put(" is "); Put(fg((i-1)*ny+ii));
    -- New_Line;

    ii := ny;
    y := yg(ii);
    kg((i-1)*ny+ii,(i-1)*ny+ii) := 1.0;
    fg((i-1)*ny+ii) := uana(x,y);
    -- Put("boundary i="&Integer'Image(i));
    -- Put(", x="); Put(x,2,3,0);
    -- Put(", ii="&Integer'Image(ii));
    -- Put(", y="); Put(y,2,3,0);
    -- Put(" is "); Put(fg((i-1)*ny+ii));
    -- New_Line;
  end loop; -- nx
  for ii in 1..ny loop
    y := yg(ii);
    i := 1;
    x := xg(i);
    kg((i-1)*ny+ii,(i-1)*ny+ii) := 1.0;
    fg((i-1)*ny+ii) := uana(x,y);
    -- Put("boundary i="&Integer'Image(i));
    -- Put(", x="); Put(x,2,3,0);
    -- Put(", ii="&Integer'Image(ii));
    -- Put(", y="); Put(y,2,3,0);
    -- Put(" is "); Put(fg((i-1)*ny+ii));
    -- New_Line;

    i := nx;
    x := xg(i);
    kg((i-1)*ny+ii,(i-1)*ny+ii) := 1.0;
    fg((i-1)*ny+ii) := uana(x,y);
    -- Put("boundary i="&Integer'Image(i));
    -- Put(", x="); Put(x,2,3,0);
    -- Put(", ii="&Integer'Image(ii));
    -- Put(", y="); Put(y,2,3,0);
    -- Put(" is "); Put(fg((i-1)*ny+ii));
    -- New_Line;
  end loop; -- ny

  -- integrate over complete domain
  Put_Line("calling gaulegf npx="&Integer'Image(npx));
  gaulegf(xmin, xmax, xx, wx, npx); -- xx and wx set up for integrations
  Put_Line("calling gaulegf npy="&Integer'Image(npy));
  gaulegf(ymin, ymax, yy, wy, npy); -- yy and wy set up for integrations
  Put_Line("xx(1)="&Real'Image(xx(1))&", xx(2)="&Real'Image(xx(2))&
           ", wx(1)="&Real'Image(wx(1))&", wx(2)="&Real'Image(wx(2)));
  Put_Line("yy(1)="&Real'Image(yy(1))&", yy(2)="&Real'Image(yy(2))&
           ", wy(1)="&Real'Image(wy(1))&", wy(2)="&Real'Image(wy(2)));
  Val := galk(xx(2),yy(2),2,2,2,2); -- four 1's in C
  Put_Line("galk(xx(2),yy(2))="&Real'Image(Val));
  Val := galf(xx(2),yy(2),2,2); -- four 1's in C
  Put_Line("galf(xx(2),yy(2))="&Real'Image(Val));

  -- compute entries in stiffness matrix
  Put_Line("compute stiffness matrix");
  for i in 2..nx-1 loop
    for j in 1..nx loop
      for ii in 2..ny-1 loop
        for jj in 1..ny loop
          -- compute integral for k(i,j)
          val := 0.0;
          for px in 1..npx loop
            for py in 1..npy loop
              val := val + wx(px)*wy(py)*galk(xx(px),yy(py),i,j,ii,jj);
            end loop;
          end loop;
          -- Put_Line("Legendre integration="&Real'Image(val)&
          --          ", at i="&Integer'Image(i)&", j="&Integer'Image(j)&
          --          ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj));
          kg((i-1)*ny+ii,(j-1)*ny+jj) := val;
        end loop; -- jj
      end loop; -- ii
    end loop; -- j
  end loop; -- i
  -- compute integral for f(i)
  for i in 2..nx-1 loop
    for ii in 2..ny-1 loop
      val := 0.0;
      for px in 1..npx loop
        for py in 1..npy loop
          val := val + wx(px)*wy(py)*galf(xx(px),yy(py),i,ii);
        end loop;
      end loop;
      -- Put_Line("Legendre integration="&Real'Image(val)&
      --          ", i="&Integer'Image(i)&", ii="&Integer'Image(ii));
      fg((i-1)*ny+ii) := val;
    end loop; -- ii
  end loop; -- i

  Put_Line("k computed stiffness matrix, see above");
  Put_Line("f computed forcing function, see above ");

  ug := simeq(kg, fg);

  avgerr := 0.0;
  maxerr := 0.0;
  Put_Line("ug computed Galerkin, Ua analytic, error ");
  for i in 1..nx loop
    for ii in 1..ny loop
      err := abs(ug((i-1)*ny+ii)-Ua((i-1)*ny+ii));
      if err>maxerr then  maxerr := err; end if;
      avgerr := avgerr + err;
      Put("ug("&Integer'Image(i)&","&Integer'Image(ii)&")=");
      Put(ug((i-1)*ny+ii),3,5,0);
      Put(", Ua="); Put(Ua((i-1)*ny+ii),3,5,0);
      Put(", err="); Put(ug((i-1)*ny+ii)-Ua((i-1)*ny+ii),3,5,2);
      New_Line;
    end loop;
  end loop;
  Write_Soln2(Ug, "pde2sin4_la.dat");
  time_now := calendar.seconds(calendar.clock);
  Put_Line("ran "&Duration'Image(Time_Now-Time_Start)&" seconds");
  Put_Line("                   maxerr="&
           Real'Image(maxerr)&", avgerr="&
           Real'image(avgerr/Real(nx*ny)));
  Put_Line("end pde2sin_la.adb");
end pde2sin_la;