-- make_abcnu_ucd.adb  generate abcnu_ucd.inp file with Dirchilet
--                   boundary values

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_Abcnu_Ucd is
  File_Name : String := "abcnu_ucd.inp";

  -- Define the region for the solution:
  -- xmin, xmin*xstep, ...
  Xmin : Real :=  -1.0;
  Ymin : Real :=  -1.0;
  Xmax : Real :=  1.0;
  Ymax : Real :=  1.0;
  Xstep : Real := 1.37;
  Ystep : Real := 1.36;
  Nx : constant := 7;
  Ny : constant := 7;
  Xinc : Real := 0.5*(Xmax-Xmin)/Real(Nx);
  Yinc : Real := 0.5*(Ymax-Ymin)/Real(Ny);
  Xg : Real_Vector(1..Nx);
  Yg : Real_Vector(1..Ny);
  Ub : Real_Vector(1..2*Nx+2*Ny);
  Xb : Real_Vector(1..2*Nx+2*Ny);
  Yb : Real_Vector(1..2*Nx+2*Ny);
  Zb : Real_Vector(1..2*Nx+2*Ny);
  X, Y : Real;
  K : 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) return Real is
  begin
    return x*x*x + 2.0*y*y*y + 3.0*x*x*y + 4.0*x*y*y + 5.0*x*y +
      6.0*x + 7.0*y + 8.0;
  end U;

begin
  Put_Line("make_abcnu_ucd.adb creating "&File_Name);
  Put_Line("nx="&Integer'Image(nx)&
           ", ny="&Integer'Image(ny));
  Put_Line("xmin="&Real'Image(Xmin)&
           ", xmax="&Real'Image(Xmax));
  Put_Line("ymin="&Real'Image(Ymin)&
           ", ymax="&Real'Image(Ymax));
  Put_Line("xstep="&Real'Image(Xstep)&
           ", ystep="&Real'Image(Ystep));
  Put_Line("xinc="&Real'Image(Xinc)&
           ", yinc="&Real'Image(Yinc));
  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 J in 1..Ny loop
     Put_Line("Y="&Real'Image(Y));
     Yg(J) := Ymin;
     Y := Y+Yinc;
     Yinc := Yinc*Ystep;
  end loop;

  K := 1;
  for I in 1..Nx loop
     Ub(K) := U(Xg(I),Yg(1));
     Xb(K) := Xg(I);
     Yb(K) := Yg(1);
     Zb(K) := 0.0;
     K := K+1;
  end loop;
  for I in 1..Nx loop
     Ub(K) := U(Xg(I),Yg(Ny));
     Xb(K) := Xg(I);
     Yb(K) := Yg(Ny);
     Zb(K) := 0.0;
     K := K+1;
  end loop;
  for J in 1..Ny loop
     Ub(K) := U(Xg(1),Yg(J));
     Xb(K) := Xg(1);
     Yb(K) := Yg(J);
     Zb(K) := 0.0;
     K := K+1;
  end loop;
  for J in 1..Ny loop
     Ub(K) := U(Xg(Nx),Yg(J));
     Xb(K) := Xg(Nx);
     Yb(K) := Yg(J);
     Zb(K) := 0.0;
     K := K+1;
  end loop;
  K := K-1;

  Create(Output, Out_File, File_Name) ;
  -- sizes n_vertices n_cells n_ndata n_cdata n_mdata
  Put(Output,K,4);
  Put(Output,K,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 line i i+1
  for I in 1..K-1 loop
     Put(Output,I,4);
     Put(Output,0,4);
     Put(Output," line ");
     Put(Output,I,4);
     Put(Output,I+1,4);
     New_Line(Output);
  end loop;
  Put(Output,K,4);
  Put(Output,0,4);
  Put(Output," line ");
  Put(Output,K,4);
  Put(Output,1,4);
  New_Line(Output);
  -- 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_abcnu_ucd.adb finished");
end Make_Abcnu_Ucd;