-- pde44e_nuderiv4d.adb
--                    forth order pde with fourth degree solution
--                    sufficient points must be available for each order,degree
--
-- given a PDE, a scatering of points in 4D, boundary values,
--
-- solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
--       utttt(x,y,z,t) + u(x,y,z,t) =  F(x,y,z,t)
-- F(x,y,z,t) = 2*(exp(x)+exp(y)+exp(z)+exp(t))
--
-- boundary conditions computed using u(x,y,z,t)
-- analytic solution is uana(x,y,z,t) = exp(x)+exp(y)+exp(z)+exp(t)
--                                   soultion, generally not known
--                                   used to compute boundary values
--                                   used for error checking
-- U(xf,yf,zf,tf) unknown at  nfree  vertices, DOF
-- U(xb,yb,zb,tb) known at  nbound  boundary vertices, nvert  total vertices
-- nuderiv4d determines cxx etc. at a subset of vertices
--
--  Uxx(xf,yf,zf,tf)= cxx(xf,yf,zf,tf)*U(xf,yf,zf,tf)+  unknown U m=70 points
--                    cxx(xfi,yfi,zfi,tfi)*U(xfi,yfi,zfi.tfi)+...+
--                             unknown U(xfi,yfi,zfi,tfi) for each xf,yf,zf,tf
--                    cxx(xbj,ybj,zbj,tbj)*U(xbj,ybj,zbj,tbj)+...
--                                                        known value
--  ...
--  f(xf,yf,zf,tf)    right hand side, known, computable
--
-- substituting in the PDE, a set of linear simultaneous equations for U(xf,yf)
--  | a11  a12  ...  a1n |   |U(xf1,yf1,zf1,tf1)| = rhs(1)  n = nfree
--  | a21  a22  ...  a2n | * |U(xf2,yf2,zf2,tf2)| = rhs(2)
--              ...
--  | an1  an2  ...  ann |   |U(xfn,yfn,zfn,tfn)| = rhs(n)
--
--
-- rhs(1) = f(xf1,yf1,zf1,tf1) - cxx(xbj,ybj,zbj,tbj)*U(xbj,ybj,zbj,tbj) -
--
--                                         ip(1)=j  any bound xbj,ybj,zbj,tbj
--
-- same code for a21, a22, ... new N.nu4d4 on next point
-- then rest of DOF points
-- solve for U(xfi,yfi)  i=1..nfree

with Ada.Text_IO;
use  Ada.Text_IO;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;
with real_Arrays;
use real_Arrays;
with simeq;
with nuderiv4d;
use nuderiv4d;

procedure pde44e_nuderiv4d is

  debug : integer := 1;
  order : integer := 4;
  nfree : integer := 0;        -- degrees of freedom-1, must be first vertices
  nvert : integer := 69;       -- number of unique vertices-1
  x : real_vector(0..nvert);
  y : real_vector(0..nvert);
  z : real_vector(0..nvert);
  t : real_vector(0..nvert);
  U : real_vector(0..nvert);          -- known boundary
  A : real_matrix(0..nfree,0..nfree); -- matrix
  rhs : real_vector(0..nfree);        -- right hand side, forcing function
  ug : real_vector(0..nfree);         -- solution computed

  nc : integer := 69;   -- present limit of nuderiv4d
  cxxxx : real_vector(0..nc);
  cyyyy : real_vector(0..nc);
  czzzz : real_vector(0..nc);
  ctttt : real_vector(0..nc);
  error, Maxerror, ap : real;
  nbound : integer := nvert-nfree;      -- number of boundary vertices
  t_start, t_init, t_kf, t_simeq, t_End : real;
  m : integer := 0;

  -- utility functions
  function max(a:real; b:real) return real is
  begin
    if a>b then
      return a;
    end if;
    return b;
  end max;

  -- test functions
  procedure gen_data4(which : integer) is
    rx : real := 0.25;
    ry : real := 0.28;
    rz : real := 0.33;
    rt : real := 0.37;
    drx : real := 0.0096;
    dry : real := 0.0094;
    drz : real := 0.0091;
    drt : real := 0.0089;
    thx : real := 0.2;
    thy : real := 1.7;
    thz : real := 2.3;
    tht : real := 3.32;
    dthx : real := 0.61;
    dthy : real := 0.72;
    dthz : real := 0.49;
    dtht : real := 0.53;
    d5th : real := 0.628318531;
    d3rd : real := 1.047197551;
    th5th : real := 1.256637061;
    r : real  := 0.2;
    theta1 : real := 0.05;
    phi1 : real := 0.1;
    psi1 : real := 0.15;
    nn : integer  := 0;
    theta, phi, psi : real  ;
  begin
    if which=0 then
      for i in 0..nvert loop
        x(i) := rx*sin(thx);
        y(i) := ry*cos(thy);
        z(i) := rz*sin(thz);
        t(i) := rt*cos(tht);
        rx := rx + drx;
        ry := ry + dry;
        rz := rz + drz;
        rt := rt + drt;
        thx := thx + dthx;
        thy := thy + dthy;
        thz := thz + dthz;
        tht := tht + dtht;
      end loop;
    else
      nn := 0;
      psi := psi1;
      for ipsi in 0..2 loop
        phi := phi1;
        for iphi in 0..4 loop
          theta := theta1;
          for itheta in 0..4 loop
            x(nn) := r*sin(psi)*sin(phi)*cos(theta);
            y(nn) := r*sin(psi)*sin(phi)*sin(theta);
            z(nn) := r*sin(psi)*cos(phi);
            t(nn) := r*cos(psi);
            nn := nn +1;
            if nn=nvert then
               return;
            end if;
            r := r + 0.011;
            theta := theta + th5th;
          end loop;
          theta1 := theta1 + 0.008;
          phi := phi + d5th;
        end loop;
        phi1 := phi1 + 0.009;
        psi := psi + d3rd;
      end loop;
    end if;
  end Gen_Data4;

  -- PDE definition functions
  function f(x: real; y: real; z: real; t: real) return real is
  begin
    return 2.0*(exp(x)+exp(y)+exp(z)+exp(t));
  end f;

  -- boundary and test function
  function uana(x : real; y : real; z : real; t : real) return real is
  begin
    return exp(x)+exp(y)+exp(z)+exp(t);
  end uana;

  -- PDE coefficients
  function aii(axxxx : real; ayyyy : real; azzzz : real; atttt : real)
               return real is
  begin
    -- solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
    --       utttt(x,y,z,t) + u(x,y,z,t) =
    return axxxx + ayyyy + azzzz + atttt + 1.0;
  end aii;

  -- PDE derivatives
  function aij(axxxx : real; ayyyy : real; azzzz : real; atttt : real)
               return real is
  begin
    -- solve uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +
    --       utttt(x,y,z,t) +    no U term
    return    axxxx +  ayyyy + azzzz + atttt;
  end aij;


begin
    Put_line("test_nuderiv4d.adb running");

    -- t_start =

    put_line("pde44e_nuderiv4d.adb running");
    put_line("  uxxxx(x,y,z,t) + uyyyy(x,y,z,t) + uzzzz(x,y,z,t) +");
    put_line("  utttt(x,y,z,t) + u(x,y,z,t) = F(x,y,z,t)");
    put_line("F(x,y,z,t) = 2*(exp(x)+exp(y)+exp(z)+exp(t))");
    put_line("uana(x,y)=exp(x)+exp(y)+exp(z)+exp(t)");
    put_line("                   solution");

    -- read in or generate data points, DOF (free) first, then boundary
    gen_data4(1);

    -- initialize vertices
    put_line("vertex, x, y, z, t");
    for i in 0..nvert loop
      U(i) := 0.0;  -- known at boundaries
      put_line(integer'image(i)&
                 " x="&real'image(x(i))&
                 ", y="&real'image(y(i)));
      put_line("    z="&real'image(z(i))&
                 ", t="&real'image(t(i)));
    end loop;

    -- clear A and rhs
    for i in 0..nfree loop
      rhs(i) := 0.0;
      for k in 0..Nfree loop
         A(i,k) := 0.0;
      end loop;
    end loop;

    -- set boundary values
    for I in 0..nvert loop
      if i>nfree then
        U(i) := uana(x(i),y(i),z(i),t(i));
        put_line("bound "&integer'image(i)&"  U(i)="&real'image(U(i)));
      end if;
    end loop;

    maxerror := 0.0;

    for point in 0..nfree loop -- for DOF
      m := nvert;
      nu4dxxxx(order, nvert, point, x, y, z, t, cxxxx);
      nu4dyyyy(order, nvert, point, x, y, z, t, cyyyy);
      nu4dzzzz(order, nvert, point, x, y, z, t, czzzz);
      nu4dtttt(order, nvert, point, x, y, z, t, ctttt);

      put_line("order="&integer'image(order)&
               ", m="&integer'image(m)&
               ", point="&integer'image(point));
      -- build row for point
      rhs(point) := f(x(point),y(point),z(point),t(point));
      A(Point,point) := A(Point,point) +
                         aii(cxxxx(point), cyyyy(point),
                             czzzz(point), ctttt(point));
      for i in 0..m loop
        if i /= point then
          ap := aij(cxxxx(i), cyyyy(i), czzzz(i), ctttt(i));

          if i<nfree then
            A(Point,i) := A(Point,i) + ap;
            if debug>0 then
              put_line("point="&integer'image(Point)&
                       ", i="&integer'image(I)&
                       ", ap="&real'image(ap));
            end if;
          else -- boundary
            rhs(point) := rhs(point) - ap*U(i);
            if debug>0 then
              put_line("point="&integer'image(Point)&
                         ", i="&integer'image(I)&
                         ", ap*U="&real'image(-ap*U(i)));
            end if;
          end if;
        end if;
      end loop; -- end i across cxxxx, cyyyy,  ...
    end loop; -- end point
    if debug>0 then
      put_line("A matrix and rhs");
      for i in 0..nfree loop
        for k in 0..nfree loop
           put_line("A("&integer'image(I)&
                    ","&integer'image(K)&
                    ")="&real'image(A(I,k)));
        end loop;
        put_line("rhs("&integer'image(I)&
                 ")="&real'image(rhs(i)));
      end loop;
      put_line(" ");
    end if;

    -- solve linear equations
    Ug := simeq(A, rhs);

    for i in 0..nfree loop
      error := abs(ug(i)-uana(x(i),y(i),z(i),t(i)));
      maxerror := max(maxerror, error);
      put_line("vertex "&integer'image(I)&
               "  x="&real'image(x(i))&
               ", y="&real'image(y(i)));
      put_line("           z="&real'image(z(i))&
               ", t="&real'image(t(i)));
      put_line("    computed="&real'image(ug(i))&
                 ", actual="&real'image(uana(x(i),y(i),z(i),t(i)))&
                 ", error = "&real'image(error));
    end loop;
    put_line("maxerror = "&Real'Image(maxerror));
    put_line(" ");
    put_line("end pde44e_nuderiv4d.adb");
end pde44e_nuderiv4d;