-- fem_bihar2dps_la.adb   biharmonic Galerkin FEM
-- solve uxxxx(x,y) + 2 uxxyy(x,y) + uyyyyy(x,y) = f(x,y)
-- f(x,y) = 4.0 * sin(x)*sin(y)
-- boundary conditions computed using ub(x,y)
-- ub(x,y) = sin(x)*sin(y)

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with 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
with rderiv;      -- function rderiv

procedure fem_bihar2dps_la is
  nx : Integer := 9;
  ny : Integer := 9;
  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);
  xmin : real := 0.0;
  xmax : real := 3.14159;
  ymin : real := 0.0;
  ymax : real := 3.14159;
  x, y, hx, hy : real;
  npx : Integer := 9;
  npy : Integer := 9;
  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);
  val, err, avgerr, maxerr : real;
  tstart, tnow, tprev : Duration;
  ib, iib : Integer;

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

  function ub(x,y:real) return real is -- analytic solution and boundary
  begin
    return sin(x)*sin(y);
  end ub;

  function galk(x,y:Real; i,j,ii,jj:integer) return real is
                                       -- Galerkin k stiffness function
  begin                                -- for this specific problem
    return (     phipppp(x,j,nx,xg)*     phi(y,jj,ny,yg) +
             2.0*  phipp(x,j,nx,xg)*   phipp(y,jj,ny,yg) +
                     phi(x,j,nx,xg)* phipppp(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;

  package Int_IO is new Integer_IO(Integer);
  use Int_IO;

  package Duration_IO is new Fixed_IO(Duration);
  use Duration_IO;

  procedure Check_Soln is -- for all interior points
    x, y, err, maxerr, avgerr, rmserr : real;
    sumerr, sumsqerr, eval_deriv : Real;
    cxx   : Real_Matrix(2..nx-1,1..nx);
    cxxxx : Real_Matrix(2..nx-1,1..nx);
    cyy   : Real_Matrix(2..ny-1,1..ny);
    cyyyy : Real_Matrix(2..ny-1,1..ny);
    tempx : Real_Vector(1..nx);
    tempy : Real_Vector(1..ny);
    tcxxyy : Real_Vector(1..ny);
    cnt : integer := 0;
  begin
    for i in 2..nx-1 loop
      rderiv(2,nx,i,hx,tempx);
      for k in 1..nx loop
        cxx(i,k) := tempx(k);
      end loop;
      rderiv(4,nx,i,hx,tempx);
      for k in 1..nx loop
        cxxxx(i,k) := tempx(k);
      end loop;
    end loop; -- i
    for ii in 2..ny-1 loop
      rderiv(2,ny,ii,hy,tempy);
      for k in 1..ny loop
        cyy(ii,k) := tempy(k);
      end loop;
      rderiv(4,ny,ii,hy,tempy);
      for k in 1..ny loop
        cyyyy(ii,k) := tempy(k);
      end loop;
    end loop; -- ii

    sumerr := 0.0;
    sumsqerr := 0.0;
    maxerr := 0.0;
    for i in 2..nx-1 loop -- full equation x
      x := xmin + Real(i-1)*hx;
      for J in 2..ny-1 loop --             y
        y := ymin + Real(j-1)*hy;
        eval_deriv := 0.0; -- no PDE u term

        for k in 1..nx loop
          eval_deriv := eval_deriv + cxxxx(i,k)*ug((k-1)*ny+j); -- PDE  uxxxx term
        end loop;
        for K in 1..ny loop
          eval_deriv := eval_deriv + cyyyy(j,k)*ug((i-1)*ny+k); -- PDE uyyyy term
        end loop;
        --  + 2*uxxyy
        for jj in 1..ny loop -- walk jj in y, computing uxx
          tcxxyy(jj) := 0.0; -- numeric  uxx part of uxxyy
          for k in 1..nx loop
            tcxxyy(jj) := tcxxyy(jj) + cxx(i,k)*ug((k-1)*ny+jj);
          end loop;
        end loop;
        for k in 1..ny loop -- uyy part of uxxyy
          eval_deriv := eval_deriv + 2.0*cyy(j,k)*tcxxyy(k); -- PDE uxxyy term
        end loop;

        err := eval_deriv-FF(x,y);
        err := abs(err);
        sumerr := sumerr + err;
        sumsqerr := sumsqerr + err*err;
        if err>maxerr then maxerr := err ; end if;
        cnt := cnt+1;
      end loop; -- y
    end loop; -- x
    rmserr := sqrt(sumsqerr/Real(cnt));
    avgerr := sumerr/Real(cnt);
    Put_Line("check_soln values");
    Put("maxerr=");
    Put(Maxerr);
    Put(", rmserr=");
    Put(Rmserr);
    Put(", avgerr=");
    Put(Avgerr);
    New_Line;
  end Check_Soln;

  procedure write_soln is
    fout : File_Type;
  begin
     Create(fout, Out_File, "fem_bihar2dps_la_ada.dat");

      Put_Line("writing fem_bihar2dps_la_ada.dat");

      for  i in 1..nx loop
        for ii in 1..ny loop
           Put(fout," ");
           put(fout,Xg(i),4,6,0);
           Put(fout," ");
           Put(fout,Yg(ii),4,6,0);
           Put(fout," ");
           Put(fout,Ug((i-1)*ny+ii),4,6,0);
           New_Line(fout);
        end loop;
        New_Line(fout);
      end loop;
      Close(fout);
      Put_Line("finished writing fem_bihar2dps_la_ada.dat");
    end write_soln;

begin
  Put_Line("fem_bihar2dps_la.adb fourth order biharmonic PDE ");
  Put_Line("solve uxxxx(x,y) + 2 uxxyy(x,y) + uyyyy(x,y) = f(x,y)");
  Put_Line("boundary conditions computed using ub(x,y)");
  Put_Line("ub(x,y) = sin(x)*sin(y)");
  Put_Line("high order shape functions used");
  Put_Line("high order Gauss-Legendre integration used");
  Put((nx-2)*(ny-2));
  Put_Line("  Degrees of Freedom");
  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));

  tstart := Calendar.Seconds(Calendar.Clock);
  tprev := tstart;

  -- initialize xg array. does not have to be uniform spacing
  hx := (xmax-xmin)/real(nx-1);
  Put_Line("x grid and boundary at ymax ");
  for i in 1..nx loop
    xg(i) := xmin+Real(i-1)*hx;
    Put("i="&Integer'Image(i)&", Ua(");
    Put(xg(i),2,5,0);
    Put(",");
    Put(ymax,2,5,0);
    Put(")=");
    Put(ub(xg(i),ymax),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 boundary at xmax");
  for ii in 1..ny loop
    yg(ii) := ymin+Real(ii-1)*hy;
    Put("ii="&Integer'Image(ii)&", Ua(");
    Put(xmax,2,5,0);
    Put(",");
    Put(yg(ii),2,5,0);
    Put(")=");
    Put(ub(xmax,yg(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) := ub(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);
    iib := 1;
    y := yg(iib);
    kg((i-1)*ny+iib,(i-1)*ny+iib) := 1.0;
    fg((i-1)*ny+iib) := ub(x,y);
    -- Put("boundary i="&Integer'Image(i));
    -- Put(", x="); Put(x,2,3,0);
    -- Put(", iib="&Integer'Image(iib));
    -- Put(", y="); Put(y,2,3,0);
    -- Put(" is "); Put(fg((i-1)*ny+iib));
    -- New_Line;

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

    ib := nx;
    x := xg(ib);
    kg((ib-1)*ny+ii,(ib-1)*ny+ii) := 1.0;
    fg((ib-1)*ny+ii) := ub(x,y);
    -- Put("boundary ib="&Integer'Image(ib));
    -- Put(", x="); Put(x,2,3,0);
    -- Put(", ii="&Integer'Image(ii));
    -- Put(", y="); Put(y,2,3,0);
    -- Put(" is "); Put(fg((ib-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); -- two 1's in C
  Put_Line("galf(xx(2),yy(2))="&Real'Image(Val));

  tnow := Calendar.Seconds(Calendar.Clock);
  Put(tnow-tprev);
  Put_Line("= wall time in seconds to initialize");
  tprev := tnow;

  -- 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");
  Put_Line("f computed, forcing function");
  tnow := Calendar.Seconds(Calendar.Clock);
  Put(tnow-tprev);
  Put_Line("= wall time in seconds to build k and f");
  tprev := tnow;

  ug := simeq(kg, fg);

  tnow := Calendar.Seconds(Calendar.Clock);
  Put(tnow-tprev);
  Put_Line("= wall time in seconds to solve k u = f");
  tprev := tnow;

  write_soln;
  check_soln;

  tnow := Calendar.Seconds(Calendar.Clock);
  Put(tnow-tprev);
  Put_Line("= wall time in seconds to write and check solution");

  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,0);
      New_Line;
    end loop;
  end loop;
  Put_Line("                   maxerr="&
           Real'Image(maxerr)&", avgerr="&
           Real'image(avgerr/Real(nx*ny)));

  tnow := Calendar.Seconds(Calendar.Clock);
  Put(tnow-tstart);
  Put_Line("= total wall time in seconds");
  Put_Line("end fem_bihar2dps_la.adb");
end fem_bihar2dps_la;