-- fem_check4th_la.adb  Galerkin FEM
-- solve  x^3 u'''' + x^2 u''' + u'' + u' + u = x^4 + 64 x^3 - x^2 - 4x + 1
-- u(0)=1, u(1)=1  0 < x < 1
-- analytic solution u(x) = x^4 - x^2 + 1
-- Gauss-Legendre integration
-- k * U = f, solve simultaneous equations for U

with Ada.Text_IO;
use  Ada.Text_IO;
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_check4th_la is
  n  : integer := 10;
  k  : Real_Matrix(1..n,1..n); -- (i,j)  Gauss-Legendre integration
  xg : Real_Vector(1..n);
  f  : Real_Vector(1..n);
  u  : Real_Vector(1..n);
  Ua : Real_Vector(1..n);
  xmin : real := 0.0;
  xmax : real := 1.0;
  h : real;
  np : integer := 24;
  xx : Real_Vector(1..np); -- for Gauss-Legendre
  ww : Real_Vector(1..np);
  val, err, avgerr, maxerr : real;

  function FF(x:real) return real is -- forcing function
  begin
    return x**4 + 64.0*x**3 - x**2 - 4.0*x + 1.0;
  end FF;

  function uana(x:real) return real is -- analytic solution and boundary
  begin
    return x**4 - x**2 + 1.0;
  end uana;

  function galk(x:Real; i,j,n:integer) return real is
                                       -- Galerkin k stiffness function
  begin                                -- for this specific problem
    return (x**3 *phipppp(x,j,n,xg) +
            x**2 * phippp(x,j,n,xg) +
            x*      phipp(x,j,n,xg) +
                     phip(x,j,n,xg) +
                      phi(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_check4th_la.c running ");
  Put_Line("Given x^3 u'''' + x^2 u''' + u'' + u' + u = x^4 + 64 x^3 - x^2 - 4x + 1");
  Put_Line("0 <= x <= 1  Boundary u(0) = 1, u(1) = 1 ");
  Put_Line("Analytic solution u(x) = x^4 - x^2 + 1 ");

  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;
    end loop;
    f(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) := 1.0

  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;
    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;
  end loop; -- i
  New_Line;

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

  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("end fem_check4th_la.adb");

end fem_check4th_la;