-- fem_checke_la.adb  Galerkin FEM
--  solve  u'(x)=e^x  u(0)=1, u(1)=e  0 < x < 1
-- initial try 10 points, then 20
-- investigate gauss-legendre, adaptive, and trapezoidal integration
-- solve for Uj
-- sum j=0,9 int x=0,1 phip(x,j)*phi(x,i) dx = int x=0,1 f(x)*phi(x,i)dx
-- k(i,j) = int x=0,1 phip(x,j)*phi(x,i) dx  (not i==0 or i==9, boundary)
-- f(i)   = int x=0,1 exp(x)*phi(x,i) dx
-- 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 aquade;       -- function aquad3 specifically for this code
use aquade;
with simeq;       -- function simeq
with gaulegf;     -- function gaulegf
with galke;       -- external galk
with galfe;       -- external galf

procedure fem_checke_la is
  k  : Real_Matrix(1..10,1..10); -- (i,j)  Gauss-Legendre integration
  ks : Real_Matrix(1..10,1..10); -- simple trapazoidal integration
  kd : Real_Matrix(1..10,1..10); -- adaptive integration
  xg : Real_Vector(1..10);
  f  : Real_Vector(1..10);
  u  : Real_Vector(1..10);
  fs : Real_Vector(1..10);
  us : Real_Vector(1..10);
  fd : Real_Vector(1..10);
  Ud : Real_Vector(1..10);
  ua : Real_Vector(1..10);
  n : integer;
  xmin : real := 0.0;
  xmax : real := 1.0;
  h, a, b : real;
  xx : Real_Vector(1..100); -- for Gauss-Legendre
  ww : Real_Vector(1..100);
  val, err, avgerr, maxerr : real;
  eps : real := 1.0e-6;
  np : integer;

  function FF(x:real) return real is -- forcing function
  begin
    return exp(x);
  end FF;

  function uana(x:real) return real is -- analytic solution and boundary
  begin
    return exp(x);
  end uana;

  -- the following two functions must be external for aquad3 using 'access
  function galk(x:Real; i,j,n:integer) return real is
                                       -- Galerkin k stiffness function
  begin                                -- for this specific problem
    return phip(x,j,n,xg)*phi(x,i,n,xg);
  end galk;

  function galf(x:Real; i:integer) return real is -- Galerkin f function
  begin
    return FF(x)*phi(x,i,n,xg);
  end galf;

  package Real_IO is new Float_IO(Real);
  use Real_IO;
begin
  Put_Line("fem_checke_la.c running ");
  Put_Line("Given du/dx = exp(x)  0<x<1  Boundary u(0) = 1, u(1) = e ");
  Put_Line("Analytic solution u(x) = exp(x) ");

  Put_Line("Try 10 points ");
  n := 10;
  h := (xmax-xmin)/real(n-1);
  Put_Line("x grid and analytic solution ");
  for i in 1..n loop
    xg(i) := xmin+Real(i-1)*h;
    Ua(i) := uana(xg(i));
    Put("i="&Integer'Image(i)&", Ua(");
    Put(xg(i),3,5,0);
    Put(")=");
    Put(Ua(i),3,5,0);
    New_Line;
  end loop;
  New_Line;

  -- initialize k and f
  for i in 1..n loop
    for j in 1..n loop
      k(i,j) := 0.0;
      ks(i,j) := 0.0;
      kd(i,j) := 0.0;
    end loop;
    f(i) := 0.0;
    fs(i) := 0.0;
    fd(i) := 0.0;
  end loop;
  -- set boundary conditions
  k(1,1) := 1.0;
  f(1) := uana(xmin);        -- u(x:=0.0) := 1.0
  k(n,n) := 1.0;
  f(n) := uana(xmax);        -- u(x:=1.0) := exp(1)
  ks(1,1) := 1.0;
  fs(1) := uana(xmin);       -- u(x:=0.0) := 1.0
  ks(n,n) := 1.0;
  fs(n) := uana(xmax);       -- u(x:=1.0) := exp(1)
  kd(1,1) := 1.0;
  fd(1) := uana(xmin);       -- u(x:=0.0) := 1.0
  kd(n,n) := 1.0;
  fd(n) := uana(xmax);       -- u(x:=1.0) := exp(1)

  a := xmin; -- integrate over complete domain
  b := xmax;
  np := 98;
  Put_Line("calling gaulegf np:="&Integer'Image(np));
  gaulegf(xmin, xmax, xx, ww, np); -- xx and ww set up for integrations
  h := (xmax-xmin)/Real(np); -- for trapezoidal integration

  -- compute entries in stiffness matrix
  Put_Line("compute stiffness matrix");
  for i in 2..n-1 loop
    for j in 1..n loop
      -- compute integral for k(i,j)
      val := 0.0;
      for p in 1..np loop
        val := val + ww(p)*galk(xx(p),i,j,n);
      end loop;
      Put_Line("Legendre integration="&Real'Image(val)&
               ", at i="&Integer'Image(i)&", j="&Integer'Image(j));
      k(i,j) := val;

      val := aquad3(galke'access, xmin, xmax, eps, i, j, n, xg);
      Put_Line("adaptive integration="&Real'Image(val)&
               ", at i="&Integer'Image(i)&", j="&Integer'Image(j));
      kd(i,j) := val;

      val := 0.0;
      for p in 1..np-1 loop
        val := val + galk(xmin+Real(p)*h,i,j,n)*h;
      end loop;
      val := val + (galk(xmin,i,j,n)+ galk(xmax,i,j,n))*h/2.0;
      Put_Line("simple   integration="&Real'Image(val)&
               ", at i="&Integer'Image(i)&", j="&Integer'Image(j));
      ks(i,j) := val;
    end loop; -- j

    -- compute integral for f(i)
    val := 0.0;
    for p in 1..np loop
      val := val + ww(p)*galf(xx(p),i);
    end loop;
    Put_Line("Legendre integration="&Real'Image(val)&", i="&Integer'Image(i));
    f(i) := val;

    val := aquad3(galfe'access, xmin, xmax, eps, i, 0, n, xg);
    Put_Line("adaptive integration="&Real'Image(val)&", i="&Integer'Image(i));
    fd(i) := val;

    val := 0.0;
    for p in 1..Np-1 loop
      val := val + galf(xmin+Real(p)*H,i)*h;
    end loop;
    val := val + (galf(xmin,i) + galf(xmax,i))*h/2.0;
    Put_Line("simple   integration="&Real'Image(val)&", i="&Integer'Image(i));
    fs(i) := val;
  end loop; -- i

  Put_Line("k computed stiffness matrix");
  for i in 1..n loop
    for j in 1..n loop
      Put_Line("i="&Integer'Image(I)&", j="&Integer'Image(J)&", k(i,j)="&
                Real'Image(k(i,j)));
    end loop;
  end loop;
  New_Line;
  Put_Line("f computed forcing function ");
  for i in 1..n loop
    Put_Line("f("&Integer'Image(i)&")="&Real'Image(f(i)));
  end loop;
  New_Line;
  u := simeq(k, f);

  avgerr := 0.0;
  maxerr := 0.0;
  Put_Line("u computed Galerkin, Ua analytic, error ");
  for i in 1..n loop
    err := abs(u(i)-Ua(i));
    if err>maxerr then  maxerr := err; end if;
    avgerr := avgerr + err;
    Put("u("&Integer'Image(I)&")=");
    Put(u(i),3,5,0);
    Put(", Ua="); Put(Ua(i),3,5,0);
    Put(", err="); Put(u(i)-Ua(i),3,5,0);
    New_Line;
  end loop;
  Put_Line("                   maxerr="&
           Real'Image(maxerr)&", avgerr="&
           Real'image(avgerr/Real(n)));
  New_Line;



  Put_Line("kd computed adaptive stiffness matrix ");
  for i in 1..n loop
    for j in 1..n loop
      Put_Line("i="&Integer'Image(I)&", j="&Integer'Image(J)&", kd(i,j)="&
                Real'Image(kd(i,j)));
    end loop;
  end loop;
  New_Line;
  Put_Line("fd computed forcing function ");
  for i in 1..n loop
    Put_Line("fd("&Integer'Image(i)&")="&Real'Image(fd(i)));
  end loop;
  New_Line;
  ud := simeq(kd, fd);

  avgerr := 0.0;
  maxerr := 0.0;
  Put_Line("ud computed Galerkin, Ua analytic, error ");
  for i in 1..n loop
    err := abs(ud(i)-Ua(i));
    if err>maxerr then  maxerr := err; end if;
    avgerr := avgerr + err;
    Put("ud("&Integer'Image(I)&")=");
    Put(ud(i),3,5,0);
    Put(", Ua="); Put(Ua(i),3,5,0);
    Put(", err="); Put(ud(i)-Ua(i),3,5,0);
    New_Line;
  end loop;
  Put_Line("                   maxerr="&
           Real'Image(maxerr)&", avgerr="&
           Real'image(avgerr/Real(n)));
  New_Line;



  Put_Line("ks computed simple stiffness matrix ");
  for i in 1..n loop
    for j in 1..n loop
      Put_Line("i="&Integer'Image(I)&", j="&Integer'Image(J)&", ks(i,j)="&
                Real'Image(ks(i,j)));
    end loop;
  end loop;
  New_Line;
  Put_Line("fs computed forcing function ");
  for i in 1..n loop
    Put_Line("fs("&Integer'Image(i)&")="&Real'Image(fs(i)));
  end loop;
  New_Line;
  us := simeq(ks, fs);

  avgerr := 0.0;
  maxerr := 0.0;
  Put_Line("us computed Galerkin, Ua analytic, error ");
  for i in 1..n loop
    err := abs(us(i)-Ua(i));
    if err>maxerr then  maxerr := err; end if;
    avgerr := avgerr + err;
    Put("us("&Integer'Image(I)&")=");
    Put(us(i),3,5,0);
    Put(", Ua="); Put(Ua(i),3,5,0);
    Put(", err="); Put(us(i)-Ua(i),3,5,0);
    New_Line;
  end loop;
  Put_Line("                   maxerr="&
           Real'Image(maxerr)&", avgerr="&
           Real'image(avgerr/Real(n)));
  New_Line;

end Fem_Checke_La;