-- fem_check22_la.adb  Galerkin FEM
-- solve uxx(x,y) + 2 uxy(x,y) + 3 uyy(x,y) +
--       4 ux(x,y)  + 5 uy(x,y)  + 6 u(x,y) = c(x,y)
-- boundary conditions computed using u(x)
-- analytic solution may be given by u(x) = x^3 + y^3 + xy + 1
-- 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 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 fem_check22_la is
  nx, ny : Integer := 4;
  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 := 1.0;
  ymin : real := 0.0;
  ymax : real := 1.0;
  x, y, hx, hy : real;
  npx , npy : Integer := 12;
  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;

  function FF(x,y:real) return real is -- forcing function
  begin
    return 6.0*x*x*x+6.0*y*y*y+12.0*x*x+15.0*y*y+6.0*x*y+11.0*x+22.0*y+8.0;
  end FF;

  function uana(x,y:real) return real is -- analytic solution and boundary
  begin
    return x*x*x+y*y*y+x*y+1.0;
  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 (      phipp(x,j,nx,xg)*  phi(y,jj,ny,yg) +
             2.0 * phip(x,j,nx,xg)* phip(y,jj,ny,yg) +
             3.0 *  phi(x,j,nx,xg)*phipp(y,jj,ny,yg) +
             4.0 * phip(x,j,nx,xg)*  phi(y,jj,ny,yg) +
             5.0 *  phi(x,j,nx,xg)* phip(y,jj,ny,yg) +
             6.0 *  phi(x,j,nx,xg)*  phi(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;
begin
  Put_Line("fem_check22_la.adb running ");
  Put_Line("solve uxx(x,y) + 2 uxy(x,y) + 3 uyy(x,y) +");
  Put_Line("      4 ux(x,y)  + 5 uy(x,y)  + 6 u(x,y) = c(x,y)");
  Put_Line("boundary conditions computed using u(x)");
  Put_Line("analytic solution may be given by u(x) = x^3 + y^3 + xy + 1");
  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));

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