-- fem_check_abc_la.adb  Galerkin FEM
-- solve a1(x,y) uxx(x,y) + b1(x,y) uxy(x,y) + c1(x,y) uyy(x,y) +
--       d1(x,y ux(x,y)   + e1(x,y) uy(x,y)  + f1(x,y) u(x,y) = c(x,y)
-- a1,b1,c1,d1,e1,f1,c and u for boundary are provided by user in abc_la.adb
-- ny, ny, xmin, xmax, ymin, ymax are provided by user in abc_la.ads
-- boundary conditions computed using u(x)
-- analytic solution may be given by u(x)
-- Gauss-Legendre integration used
-- solve for Uij numerical approximation U(x_i,y_j) of u(x_i,y_j)
-- kg * U = fg, solve simultaneous equations for U , ug.

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
with abc_la;      -- problem definition
use  abc_la;

procedure fem_check_abc_la is
  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);
  x, y, hx, hy : real;
  xx : Real_Vector(1..100); -- for Gauss-Legendre
  wx : Real_Vector(1..100);
  yy : Real_Vector(1..100); -- for Gauss-Legendre
  wy : Real_Vector(1..100);
  npx , npy, i, ii : Integer;
  val, err, avgerr, maxerr : real;

  function galk(x,y:Real; i,j,ii,jj:integer) return real is
                                       -- Galerkin k stiffness function
  begin                                -- for this specific problem
    return (a1(x,y)*phipp(x,j,nx,xg)*  phi(y,jj,ny,yg) +
            b1(x,y)* phip(x,j,nx,xg)* phip(y,jj,ny,yg) +
            c1(x,y)*  phi(x,j,nx,xg)*phipp(y,jj,ny,yg) +
            d1(x,y)* phip(x,j,nx,xg)*  phi(y,jj,ny,yg) +
            e1(x,y)*  phi(x,j,nx,xg)* phip(y,jj,ny,yg) +
            f1(x,y)*  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 c(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_check_abc_la.adb running ");
  Put_Line("solve a1(x,y) uxx(x,y) + b1(x,y) uxy(x,y) + c1(x,y) uyy(x,y) +");
  Put_Line("      d1(x,y) ux(x,y)  + e1(x,y) uy(x,y)  + f1(x,y) u(x,y) = c(x,y)");
  Put_Line("boundary conditions computed using u(x)");
  Put_Line("user provides a1,b1,c1,d1,e1,f1,c,u in abc_la.adb");
  Put_Line("user provides nx,ny,xmin,xmax,ymin,ymax in abc_la.ads");
  Put_Line("analytic solution may be given by u(x) with ifcheck>0");
  Put_Line("Gauss-Legendre integration used");
  New_Line;
  Put("xmin="); Put(xmin,3,2,0);
  Put(", xmax="); Put(xmax,3,2,0);
  Put(", ymin="); Put(ymin,3,2,0);
  Put(", ymax="); Put(ymax,3,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) := u(xg(i),ymin);
    Put("i="&Integer'Image(i)&", Ua(");
    Put(xg(i),3,5,0);
    Put(",");
    Put(ymin,3,5,0);
    Put(")=");
    Put(Ua(i),4,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) := u(xmin,yg(ii));
    Put("ii="&Integer'Image(ii)&", Ua(");
    Put(xmin,3,5,0);
    Put(",");
    Put(yg(ii),3,5,0);
    Put(")=");
    Put(Ua(ii),4,5,0);
    New_Line;
  end loop;
  New_Line;

  -- put solution in Ua vector
  if ifcheck>0 then
    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) := u(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;
  end if; -- ifcheck

  -- 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) := u(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) := u(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) := u(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) := u(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
  npx := 48;
  npy := 48;
  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

  -- 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;
          if ifcheck>1 then
            Put_Line("Legendre integration="&Real'Image(val)&
                     ", at i="&Integer'Image(i)&", j="&Integer'Image(j)&
                     ", ii="&Integer'Image(ii)&", jj="&Integer'Image(jj));
          end if; -- ifcheck
          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;
      if ifcheck>1 then
        Put_Line("Legendre integration="&Real'Image(val)&
                 ", i="&Integer'Image(i)&", ii="&Integer'Image(ii));
      end if; -- ifcheck
      fg((i-1)*ny+ii) := val;
    end loop; -- ii
  end loop; -- i

  ug := simeq(kg, fg);

  if ifcheck>0 then
    avgerr := 0.0;
    maxerr := 0.0;
    Put_Line("solution 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)));
  else -- no ifcheck
    Put_Line("solution ug computed by Galerkin method");
    for i in 1..nx loop
      for ii in 1..ny loop
        Put("ug("&Integer'Image(i)&","&Integer'Image(ii)&")=");
        Put(ug((i-1)*ny+ii),3,5,0);
        New_Line;
      end loop;
    end loop;
  end if; -- ifcheck
  Put_Line("end fem_check_abc_la.adb");
end fem_check_abc_la;