-- make_nu3_ucd.adb  generate nu3_ucd.inp file with Dirchilet
--                   boundary values
-- pde3b.adb  3D non equal h boundary value PDE
--  u(x,y,z) = x^3  + y^3  +z^3  + 2x^2y + 3xy^2 + 4z^2x + 5y + 6

with Ada.Text_IO;
use  Ada.Text_IO;
with Real_Arrays;
use  Real_Arrays;
with Ada.Numerics.Long_Elementary_Functions;
use  Ada.Numerics.Long_Elementary_Functions;

procedure Make_Nu3_Ucd is
  File_Name : String := "nu3_ucd.inp";

  -- Define the region for the solution:
  -- xmin, xmin*xstep, ...
  Xmin : Real :=  -1.0;
  Ymin : Real :=  -1.0;
  Zmin : Real :=  -1.0;
  Xmax : Real :=  1.0;
  Ymax : Real :=  1.0;
  Zmax : Real :=  1.0;
  Xstep : Real := 1.48;
  Ystep : Real := 1.37;
  Zstep : Real := 1.28;
  nx:constant   :=     6;
  ny:constant   :=     7;
  nz:constant   :=     8;
  Xinc : Real := 0.5*(Xmax-Xmin)/Real(Nx);
  Yinc : Real := 0.5*(Ymax-Ymin)/Real(Ny);
  Zinc : Real := 0.5*(Zmax-Zmin)/Real(Nz);
  Xg : Real_Vector(1..Nx);
  Yg : Real_Vector(1..Ny);
  Zg : Real_Vector(1..Nz);
  Ub : Real_Vector(1..2*Nx*Ny+2*Ny*Nz+2*Nx*Nz);
  Xb : Real_Vector(1..2*Nx*Ny+2*Ny*Nz+2*Nx*Nz);
  Yb : Real_Vector(1..2*Nx*Ny+2*Ny*Nz+2*Nx*Nz);
  Zb : Real_Vector(1..2*Nx*Ny+2*Ny*Nz+2*Nx*Nz);
  type Cell is record
     N1, N2, N3, N4 : Integer;
  end record;
  C : array(1..2*(Nx-1)*(Ny-1)+2*(Ny-1)*(Nz-1)+2*(Nx-1)*(Nz-1)) of Cell;
  X, Y, Z : Real;
  K, Kc : Integer;

  package Int_Io is new Integer_Io(Integer);
  use Int_Io;
  package Real_Io is new Float_Io(Real);
  use Real_Io;
  Output : File_Type;

  function u(x:real; y:real; z:real) return real is
           -- solution for boundary and check
  begin
    return x*x*x + y*y*y + z*z*z + 2.0*x*x*y + 3.0*x*y*y +
           4.0*z*z*x + 5.0*y +6.0;
  end u;

begin
  Put_Line("make_nu3_ucd.adb creating "&File_Name);
  Put_Line("nx="&Integer'Image(nx)&", ny="&Integer'Image(ny)&
           ", nz="&Integer'Image(nz));
  Put_Line("xmin="&Real'Image(Xmin)&", ymin="&Real'Image(Ymin)&
           ", zmin="&Real'Image(Zmin));
  Put_Line("xstep="&Real'Image(Xstep)&", ystep="&Real'Image(Ystep)&
           ", zstep="&Real'Image(Zstep));
  Put_Line("xinc="&Real'Image(Xinc)&", yinc="&Real'Image(Yinc)&
           ", zinc="&Real'Image(Zinc));
  X := Xmin;
  for I in 1..Nx loop
     Put_Line("X="&Real'Image(X));
     Xg(I) := X;
     X := X+Xinc;
     Xinc := Xinc*Xstep;
  end loop;
  Y := Ymin;
  for I in 1..Ny loop
     Put_Line("Y="&Real'Image(Y));
     Yg(I) := Y;
     Y := Y+Yinc;
     Yinc := Yinc*Ystep;
  end loop;
  Z := Zmin;
  for I in 1..Nz loop
     Put_Line("Z="&Real'Image(Z));
     Zg(I) := Z;
     Z := Z+Zinc;
     Zinc := Zinc*Zstep;
  end loop;

  K := 1; -- count vertices
  Kc := 1; -- count cells
  for I in 1..Nx loop
    for J in 1..Ny loop
       Ub(K) := U(Xg(I),Yg(J),Zg(1));
       Xb(K) := Xg(I);
       Yb(K) := Yg(J);
       Zb(K) := Zg(1);
       if I>1 and J>1 then
          C(Kc).N1 := K-Ny-1;
          C(Kc).N2 := K-Ny;
          C(Kc).N3 := K;
          C(Kc).N4 := K-1;
          Kc := Kc+1;
       end if;
       K := K+1;
    end loop;
  end loop;
  for I in 1..Nx loop
    for J in 1..Ny loop
      Ub(K) := U(Xg(I),Yg(J),Zg(Nz));
      Xb(K) := Xg(I);
      Yb(K) := Yg(J);
      Zb(K) := Zg(Nz);
       if I>1 and J>1 then
          C(Kc).N1 := K-Ny-1;
          C(Kc).N2 := K-Ny;
          C(Kc).N3 := K;
          C(Kc).N4 := K-1;
          Kc := Kc+1;
       end if;
      K := K+1;
    end loop;
  end loop;

  for I in 1..Nx loop
    for J in 1..Nz loop
      Ub(K) := U(Xg(I),Yg(1),Zg(J));
      Xb(K) := Xg(I);
      Yb(K) := Yg(1);
      Zb(K) := Zg(J);
       if I>1 and J>1 then
          C(Kc).N1 := K-Nz-1;
          C(Kc).N2 := K-Nz;
          C(Kc).N3 := K;
          C(Kc).N4 := K-1;
          Kc := Kc+1;
       end if;
      K := K+1;
    end loop;
  end loop;
  for I in 1..Nx loop
    for J in 1..Nz loop
      Ub(K) := U(Xg(I),Yg(Ny),Zg(J));
      Xb(K) := Xg(I);
      Yb(K) := Yg(Ny);
      Zb(K) := Zg(J);
       if I>1 and J>1 then
          C(Kc).N1 := K-Nz-1;
          C(Kc).N2 := K-Nz;
          C(Kc).N3 := K;
          C(Kc).N4 := K-1;
          Kc := Kc+1;
       end if;
      K := K+1;
    end loop;
  end loop;

  for I in 1..Ny loop
    for J in 1..Nz loop
      Ub(K) := U(Xg(1),Yg(I),Zg(J));
      Xb(K) := Xg(1);
      Yb(K) := Yg(I);
      Zb(K) := Zg(J);
       if I>1 and J>1 then
          C(Kc).N1 := K-Nz-1;
          C(Kc).N2 := K-Nz;
          C(Kc).N3 := K;
          C(Kc).N4 := K-1;
          Kc := Kc+1;
       end if;
      K := K+1;
    end loop;
  end loop;
  for I in 1..Ny loop
    for J in 1..Nz loop
      Ub(K) := U(Xg(Nx),Yg(I),Zg(J));
      Xb(K) := Xg(Nx);
      Yb(K) := Yg(I);
      Zb(K) := Zg(J);
       if I>1 and J>1 then
          C(Kc).N1 := K-Nz-1;
          C(Kc).N2 := K-Nz;
          C(Kc).N3 := K;
          C(Kc).N4 := K-1;
          Kc := Kc+1;
       end if;
      K := K+1;
    end loop;
  end loop;

  K := K-1;
  Kc := Kc-1;

  Create(Output, Out_File, File_Name) ;
  -- sizes n_vertices n_cells n_ndata n_cdata n_mdata
  Put(Output,K,4);
  Put(Output,Kc,4);
  Put(Output,1,4);
  Put(Output,0,4);
  Put(Output,0,4);
  New_Line(Output);
  -- vertices i x y z
  for I in 1..K loop
     Put(Output,I,4);
     Put(Output," "); Put(Output,Xb(I),4,6,0);
     Put(Output," "); Put(Output,Yb(I),4,6,0);
     Put(Output," "); Put(Output,Zb(I),4,6,0);
     New_Line(Output);
  end loop;
  -- cells i 0 quad n1 n2 n3 n4
  for I in 1..Kc loop
     Put(Output,I,4);
     Put(Output,0,4);
     Put(Output," quad ");
     Put(Output,C(I).N1,4);
     Put(Output,C(I).N2,4);
     Put(Output,C(I).N3,4);
     Put(Output,C(I).N4,4);
     New_Line(Output);
  end loop;

  -- node data counts
  Put(Output,1,4);
  Put(Output,1,4);
  New_Line(Output);
  Put_Line(Output,"dirchilet boundary values");
  for I in 1..K loop
     Put(Output,I,4);
     Put(Output," "); Put(Output,Ub(I),4,6,0);
     New_Line(Output);
  end loop;

  Close(Output);
  Put_Line("make_nu3_ucd.adb finished");
end Make_Nu3_Ucd;