-- fem_nl12_la.adb  Galerkin FEM
--
--
-- solve Dc*ux(x) + E(x)*(-ux(x)/u(x)^2 - Fc*uxx(x) = f(x)
--       f(x) = -2
--       Dc = 1
--       Fc = 1
--       E(X) = (x^2+1)^2
--       Analytic Solution U(X)=x^2+1
--       Inside Grid, Gauss Legendre Integration
--
--
-- solve for Uj
-- k(i,j) = int x=0,1 galk(i,j) dx  (not i==0 or i==9, boundary)
-- f(i)   = int x=0,1 FF(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
use  laphi;
with simeq;       -- function simeq
with gaulegf;     -- function gaulegf

procedure fem_nl12_la is
  nx : Integer := 10;
  k  : Real_Matrix(1..nx,1..nx); -- (i,j)  Gauss-Legendre integration
  xg : Real_Vector(1..nx);
  f  : Real_Vector(1..nx);
  u  : Real_Vector(1..nx);
  ua : Real_Vector(1..nx);
  n : integer;
  xmin : real := 1.0;
  xmax : real := 2.0;
  h : 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;

  -- definition of PDE to be solved
  Dc : Real := 1.0;
  Fc : Real := 1.0;
  function E(x:real) return real is -- coefficient function
  begin
    return (x*x+1.0)*(x*x+1.0);
  end E;

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

  function uana(x:real) return real is -- analytic solution and boundary
  begin
    return x*x+1.0;
  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 (Dc*phip(x,j,n,xg)
             -E(x)*phip(x,j,n,xg)/(phi(x,j,n,xg)*phi(x,j,n,xg))
             -Fc*phipp(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_nl12_la.c running ");
  Put_Line("solve Dc*ux(x) + E(x)*(-ux(x)/u(x)^2 - Fc*uxx(x) = f(x)");
  Put_Line("      f(x) = -2");
  Put_Line("      Dc = 1");
  Put_Line("      Fc = 1");
  Put_Line("      E(X) = (x^2+1)^2");
  Put_Line("Analytic solution u(x) = x^2+1 ");
  Put_Line("xmin="&Real'Image(Xmin)&", xmax="&Real'Image(Xmax));

  Put_Line("nx="&Integer'Image(nx)&" points including boundary");
  n := 10;
  h := (xmax-xmin)/real(nx-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(xmin)
  k(n,n) := 1.0;
  f(n) := uana(xmax);        -- u(xmax)

  np := 24;  -- integrate over complete domain
  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; -- p
      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

  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;

end fem_nl12_La;