-- make_abc_ucd.adb  generate abc_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_Abc_Ucd is
  File_Name : String := "abc_ucd.inp";

  -- Define the region for the solution:
  -- xmin <= x <= xmax, ymin <= y <= ymax  using  nx by ny points
  xmin:constant :=  -1.0;
  xmax:constant :=   1.0;
  ymin:constant :=  -1.0;
  ymax:constant :=   1.0;
  Nx : constant :=     7;
  Ny : constant :=     7;
  hx:Real := (xmax-xmin)/Real(nx-1); -- x step
  hy:Real := (ymax-ymin)/Real(ny-1); -- y step
  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);
  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_abc_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("hx="&Real'Image(hx)&
           ", hy="&Real'Image(hy));

  for I in 1..Nx loop
     Xg(I) := Xmin+Real(I-1)*Hx;
  end loop;
  for J in 1..Ny loop
     Yg(J) := Ymin+Real(J-1)*Hy;
  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_abc_ucd.adb finished");
end Make_Abc_Ucd;