with
  text_io,
  ada.integer_text_io,
  ada.long_float_text_io,
  ada.numerics.long_elementary_functions,
  terms,
  nlp;

procedure iconvex2 is

  subtype fmat3x1 is terms.fmat1 (1 .. 3);

  subtype fmat6x1 is terms.fmat1 (1 .. 6);

  type elem1_rec is
    record
      glob   : terms.imat1 (1 .. 2);
      u      : terms.fmat1 (1 .. 3);
      area   : long_float;
      length : long_float;
    end record;

  type hyperelastic is (neo, mooney);

  type mat_rec (h : hyperelastic) is
    record
      hyper : hyperelastic;

      case h is
        when neo =>
          mu : long_float;

        when mooney =>
          mu10 : long_float;
          mu01 : long_float;
      end case;
    end record;

  type mat_ptr is access mat_rec;

  type elem1_vec is array (integer range <>) of elem1_rec;

  type elem1_ptr is access elem1_vec;

  optimal                           : boolean;
  nn,n,maxls,minor,maxminor,maxlist : integer;
  lstol,mntol,mnref,maxnorm,f       : long_float;
  index                             : terms.imat2_ptr;
  xo,yo,zo,x,g                      : terms.fmat1_ptr;
  known,force                       : terms.fmat2_ptr;
  elem1                             : elem1_ptr;
  mat                               : mat_ptr;

  procedure scan (io : in text_io.file_type) is

    ascii : character;

  begin
    loop
      text_io.get(io,ascii);
      text_io.skip_line(io);
      exit when ascii = ':';
    end loop;
  end scan;

  procedure input is

    use
      text_io,
      ada.integer_text_io,
      ada.long_float_text_io;

    e,ne,glob,axis,nk,h : integer;
    io                  : file_type;

  begin
    open(io, in_file, "input.txt");

    -- The iterations terminate if the infinity norm of the
    -- gradient becomes less than or equal to MAXNORM, which is
    -- calculated as:

    -- MAXNORM = MNTOL * MNREF

    -- If MNTOL is specified as a non-positive value, it is
    -- replaced by 0.001.

    -- If MNREF is specified as a non-positive value, it is
    -- replaced by the infinite norm of the gradient, evaluated at
    -- the starting point.

    scan(io);
    get(io,lstol);
    get(io,maxls);
    get(io,mntol);
    get(io,mnref);
    get(io,maxminor);
    get(io,maxlist);
    skip_line(io);

    if mntol <= 0.0 then
      mntol := 0.001;
    end if;

    -- material
    scan(io);
    get(io,h);

    case h is
      when 1 =>
        mat := new mat_rec (neo);
        mat.hyper := neo;
        get(io,mat.mu);

      when 2 =>
        mat := new mat_rec (mooney);
        mat.hyper := mooney;
        get(io,mat.mu10);
        get(io,mat.mu01);

      when others =>
        null;
    end case;

    skip_line(io);

    -- nodes
    scan(io);
    get(io,nn);
    skip_line(io);

    -- initiate xo, yo, zo
    xo := new terms.fmat1 (1 .. nn);
    yo := new terms.fmat1 (1 .. nn);
    zo := new terms.fmat1 (1 .. nn);

    skip_line(io);
    for k in 1 .. nn loop
      get(io,glob);
      get(io,xo(glob));
      get(io,yo(glob));
      get(io,zo(glob));
      skip_line(io);
    end loop;

    -- elements 1
    scan(io);
    get(io,ne);
    skip_line(io);

    elem1 := new elem1_vec (1 .. ne);

    if ne > 0 then
      skip_line(io);
      for k in elem1'range loop
        get(io,e);
        get(io,elem1(e).glob(1));
        get(io,elem1(e).glob(2));
        get(io,elem1(e).area);
        skip_line(io);
      end loop;
    end if;

    -- initiate index
    index := new terms.imat2 (1 .. nn, 1 .. 3);
    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        index(glob,axis) := 0;
      end loop;
    end loop;

    -- initiate known
    known := new terms.fmat2 (1 .. nn, 1 .. 3);
    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        known(glob,axis) := 0.0;
      end loop;
    end loop;

    -- known displacements
    scan(io);
    get(io,nk);
    skip_line(io);

    if nk > 0 then
      skip_line(io);
      for k in 1 .. nk loop
        get(io,glob);
        get(io,axis);
        get(io,known(glob,axis));
        skip_line(io);
        index(glob,axis) := 1;
      end loop;
    end if;

    -- initiate force
    force := new terms.fmat2 (1 .. nn, 1 .. 3);
    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        force(glob,axis) := 0.0;
      end loop;
    end loop;

    -- nodal forces
    scan(io);
    get(io,nk);
    skip_line(io);

    if nk > 0 then
      skip_line(io);
      for k in 1 .. nk loop
        get(io,glob);
        get(io,axis);
        get(io,force(glob,axis));
        skip_line(io);
      end loop;
    end if;

    close(io);
  end input;

  procedure element1 (elem1     : in elem1_rec;
                      fe,stress : out long_float;
                      x         : in terms.fmat1_ptr;
                      ge        : out fmat6x1) is

    use
      ada.numerics.long_elementary_functions;

    glob,i                      : integer;
    psi,psi1,psi2,c1,c2,ksi,aux : long_float;
    z,v                         : terms.fmat1 (1 .. 3);
    xe                          : terms.fmat2 (1 .. 2, 1 .. 3);

  begin
    -- global displacements
    for loc in 1 .. 2 loop
      glob := elem1.glob(loc);
      for axis in 1 .. 3 loop
        i := index(glob,axis);
        if i > n then
          xe(loc,axis) := known(glob,axis);
        else
          xe(loc,axis) := x(i);
        end if;
      end loop;
    end loop;

    -- z
    for axis in 1 .. 3 loop
      z(axis) := (xe(2,axis) - xe(1,axis)) / elem1.length;
    end loop;

    -- ksi
    ksi := 0.0;
    for axis in 1 .. 3 loop
      ksi := ksi + z(axis) * (2.0 * elem1.u(axiS) + z(axis));
    end loop;

    -- Invariants of the right Cauchy-Green deformation tensor
    c1 := (1.0 + ksi) + 2.0 / sqrt(1.0 + ksi);
    c2 := (1.0 + ksi) ** 2 + 2.0 / (1.0 + ksi);

    -- strain energy per unit undeformed volume
    case mat.hyper is
      when neo =>
        psi := mat.mu * c1 / 2.0;
        psi1 := mat.mu / 2.0;
        psi2 := 0.0;

      when mooney =>
        psi := mat.mu10 * c1 + mat.mu01 * (c1 * c1 - c2) / 2.0;
        psi1 := mat.mu10 + mat.mu01 * c1;
        psi2 := -mat.mu01 / 2.0;
    end case;

    -- objective
    fe := psi * elem1.area * elem1.length;

    -- Cauchy stress
    aux := sqrt(1.0 + ksi);
    stress := 2.0 * ksi * (aux + 1.0 / (aux + 1.0)) / aux * (psi1 + 2.0 * psi2 * (1.0 + ksi + 1.0 / aux));

    -- gradient
    aux := elem1.area / (1.0 + ksi) * stress;
    for axis in 1 .. 3 loop
      v(axis) := aux * (elem1.u(axis) + z(axis));
    end loop;

    ge(1) := -v(1);
    ge(2) := -v(2);
    ge(3) := -v(3);
    ge(4) := v(1);
    ge(5) := v(2);
    ge(6) := v(3);
  end element1;

  procedure objective (n : in integer;
                       f : out long_float;
                       x : in terms.fmat1_ptr;
                       g : in terms.fmat1_ptr) is

    glob,i,s  : integer;
    fe,stress : long_float;
    ge1       : fmat6x1;

  begin
    -- initiate objective, gradient
    f := 0.0;
    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        i := index(glob,axis);
        if i <= n then
          f := f - force(glob,axis) * x(i);
          g(i) := -force(glob,axis);
        end if;
      end loop;
    end loop;

    -- elements 1
    for e in elem1'range loop
      element1(elem1(e),fe,stress,x,ge1);

      -- update objective
      f := f + fe;

      -- update gradient
      s := 0;
      for loc in 1 .. 2 loop
        glob := elem1(e).glob(loc);
        for axis in 1 .. 3 loop
          s := s + 1;
          i := index(glob,axis);
          if i <= n then
            g(i) := g(i) + ge1(s);
          end if;
        end loop;
      end loop;
    end loop;
  end objective;

  package qn is new nlp.qnewton(objective);

  procedure start is

    use
      ada.numerics.long_elementary_functions;

    glob1,glob2,s1,s2 : integer;
    norm              : long_float;

  begin
    -- number of unknowns
    n := 3 * nn;
    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        if index(glob,axis) = 1 then
          n := n - 1;
        end if;
      end loop;
    end loop;

    -- new sequence
    s1 := 0;
    s2 := n;
    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        if index(glob,axis) = 0 then
          s1 := s1 + 1;
          index(glob,axis) := s1;
        else
          s2 := s2 + 1;
          index(glob,axis) := s2;
        end if;
      end loop;
    end loop;

    -- elements 1
    for e in elem1'range loop
      glob1 := elem1(e).glob(1);
      glob2 := elem1(e).glob(2);

      -- u
      elem1(e).u(1) := xo(glob2) - xo(glob1);
      elem1(e).u(2) := yo(glob2) - yo(glob1);
      elem1(e).u(3) := zo(glob2) - zo(glob1);

      norm := sqrt(elem1(e).u(1) ** 2 + elem1(e).u(2) ** 2 + elem1(e).u(3) ** 2);

      elem1(e).u(1) := elem1(e).u(1) / norm;
      elem1(e).u(2) := elem1(e).u(2) / norm;
      elem1(e).u(3) := elem1(e).u(3) / norm;

      -- length
      elem1(e).length := norm;
    end loop;

    -- initiate x
    x := new terms.fmat1 (1 .. n);
    for k in 1 .. n loop
      x(k) := 0.0;
    end loop;

    -- initiate g
    g := new terms.fmat1 (1 .. n);
    for k in 1 .. n loop
      g(k) := 0.0;
    end loop;

    objective(n,f,x,g);

    -- maxnorm
    if mnref <= 0.0 then
      maxnorm := mntol * nlp.norm(n,g);
    else
      maxnorm := mntol * mnref;
    end if;
  end start;

  procedure vertex (glob : in integer;
                    vx   : out fmat3x1) is

    i : integer;

  begin
    vx(1) := xo(glob);
    vx(2) := yo(glob);
    vx(3) := zo(glob);

    for axis in 1 .. 3 loop
      i := index(glob,axis);
      if i <= n then
        vx(axis) := vx(axis) + x(i);
      else
        vx(axis) := vx(axis) + known(glob,axis);
      end if;
    end loop;
  end vertex;

  procedure autocad is

    use
      text_io,
      ada.integer_text_io,
      ada.long_float_text_io;

    glob      : integer;
    fe,stress : long_float;
    ge1       : fmat6x1;
    vx        : fmat3x1;
    io        : file_type;

  begin
    create(io, out_file, "autocad.scr");

    -- set running object snap mode to NONE
    put(io,"osmode 0");
    new_line(io);

    -- elements 1
    if elem1'last > 0 then
      put(io,"layer make Initial_1 ");
      new_line(io);

      for e in elem1'range loop
        put(io,"3dpoly ");
        for loc in 1 .. 2 loop
          glob := elem1(e).glob(loc);
          put(io,xo(glob),1,8,3);
          put(io,",");
          put(io,yo(glob),1,8,3);
          put(io,",");
          put(io,zo(glob),1,8,3);
          put(io," ");
        end loop;
        new_line(io);
      end loop;

      put(io,"layer make Final_1 ");
      new_line(io);

      for e in elem1'range loop
        element1(elem1(e),fe,stress,x,ge1);
        if stress = 0.0 then
          put(io,"-color");
          new_line(io);
          put(io,"white");
          new_line(io);
        elsif stress < 0.0 then
          put(io,"-color");
          new_line(io);
          put(io,"red");
          new_line(io);
        elsif stress > 0.0 then
          put(io,"-color");
          new_line(io);
          put(io,"blue");
          new_line(io);
        end if;

        put(io,"3dpoly ");
        for loc in 1 .. 2 loop
          glob := elem1(e).glob(loc);
          vertex(glob,vx);
          put(io,vx(1),1,8,3);
          put(io,",");
          put(io,vx(2),1,8,3);
          put(io,",");
          put(io,vx(3),1,8,3);
          put(io," ");
        end loop;
        new_line(io);
      end loop;

      put(io,"-color");
      new_line(io);
      put(io,"white");
      new_line(io);

      put(io,"layer off Initial_1 ");
      new_line(io);
    end if;

    put(io,"layer set 0 ");
    new_line(io);

    put(io,"zoom ext ext");
    new_line(io);

    close(io);
  end autocad;

  procedure putchar (io : in text_io.file_type;
                     c  : in character;
                     n  : in integer) is

  begin
    for k in 1 .. n loop
      text_io.put(io,c);
    end loop;
  end putchar;

  procedure output is

    use
      text_io,
      ada.integer_text_io,
      ada.long_float_text_io;

    i,glob,last,s : integer;
    fe,stress     : long_float;
    ge1           : fmat6x1;
    io            : file_type;

  begin
    create(io, out_file, "output.txt");

    putchar(io,'-',63);
    new_line(io);

    put(io,"| Minor Iterations                            | ");
    put(io,minor,13);
    put(io," |");
    new_line(io);

    put(io,"| Objective                                   | ");
    put(io,f,2,6,3);
    put(io," |");
    new_line(io);

    put(io,"| Allowable Maximum Gradient                  | ");
    put(io,maxnorm,2,6,3);
    put(io," |");
    new_line(io);

    putchar(io,'-',63);
    new_line(io);

    -- displacements
    new_line(io);
    new_line(io);
    putchar(io,'-',59);
    new_line(io);
    put(io," Node           Displ X           Displ Y           Displ Z");
    new_line(io);
    putchar(io,'-',59);
    new_line(io);

    new_line(io);

    for glob in 1 .. nn loop
      put(io,glob,5);
      for axis in 1 .. 3 loop
        putchar(io,' ',5);
        i := index(glob,axis);
        if i > n then
          put(io,known(glob,axis),2,6,3);
        else
          put(io,x(i),2,6,3);
        end if;
      end loop;
      new_line(io);
    end loop;

    -- support reactions
    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        force(glob,axis) := -force(glob,axis);
      end loop;
    end loop;

    for e in elem1'range loop
      element1(elem1(e),fe,stress,x,ge1);

      s := 0;
      for loc in 1 .. 2 loop
        glob := elem1(e).glob(loc);
        for axis in 1 .. 3 loop
          s := s + 1;
          force(glob,axis) := force(glob,axis) + ge1(s);
        end loop;
      end loop;
    end loop;

    new_line(io);
    new_line(io);
    putchar(io,'-',59);
    new_line(io);
    put(io," Node           Force X           Force Y           Force Z");
    new_line(io);
    putchar(io,'-',59);
    new_line(io);

    for glob in 1 .. nn loop
      last := 0;
      for axis in 1 .. 3 loop
        if index(glob,axis) > n and abs force(glob,axis) > maxnorm then
          last := axis;
        end if;
      end loop;

      if last /= 0 then
        put(io,glob,5);
        for axis in 1 .. last loop
          if index(glob,axis) > n and abs force(glob,axis) > maxnorm then
            putchar(io,' ',5);
            put(io,force(glob,axis),2,6,3);
          else
            putchar(io,' ',18);
          end if;
        end loop;
        new_line(io);
      end if;
    end loop;

    -- Cauchy stress
    new_line(io);
    new_line(io);
    putchar(io,'-',23);
    new_line(io);
    put(io," Elem            Stress");
    new_line(io);
    putchar(io,'-',23);
    new_line(io);

    for e in elem1'range loop
      element1(elem1(e),fe,stress,x,ge1);

      put(io,e,5);
      putchar(io,' ',5);
      put(io,stress,2,6,3);
      new_line(io);
    end loop;

    close(io);
  end output;

  procedure abaqus is

    use
      text_io,
      ada.integer_text_io,
      ada.long_float_text_io;

    nsubst : integer :=  10;
    neqit  : integer := 100;
    io     : file_type;

  begin
    create(io, out_file, "abaqusfile.inp");

    put(io,"**");
    new_line(io);
    put(io,"** NODE COORDINATES");
    new_line(io);
    put(io,"*NODE");
    new_line(io);

    for glob in 1 .. nn loop
      put(io,glob,5);
      put(io,", ");
      put(io,xo(glob),2,6,3);
      put(io,", ");
      put(io,yo(glob),2,6,3);
      put(io,", ");
      put(io,zo(glob),2,6,3);
      new_line(io);
    end loop;

    put(io,"**");
    new_line(io);
    put(io,"** ELEMENTS");
    new_line(io);

    for e in elem1'range loop
      put(io,"*ELEMENT, TYPE=T3D2, ELSET=SET");
      put(io,e,0);
      new_line(io);
      put(io,e,5);
      put(io,", ");
      put(io,elem1(e).glob(1),5);
      put(io,", ");
      put(io,elem1(e).glob(2),5);
      new_line(io);
    end loop;

    put(io,"**");
    new_line(io);
    put(io,"** MATERIAL DEFINITION");
    new_line(io);
    put(io,"*MATERIAL, NAME=MAT");
    new_line(io);

    case mat.hyper is
      when neo =>
        put(io,"*HYPERELASTIC, NEO HOOKE");
        new_line(io);
        put(io,(mat.mu / 2.0),2,6,3);
        put(io,", 0.000000E+00");
        new_line(io);

      when mooney =>
        put(io,"*HYPERELASTIC, MOONEY-RIVLIN");
        new_line(io);
        put(io,mat.mu10,2,6,3);
        put(io,", ");
        put(io,mat.mu01,2,6,3);
        put(io,", 0.000000E+00");
        new_line(io);
    end case;

    put(io,"**");
    new_line(io);
    put(io,"** ELEMENT PROPERTIES");
    new_line(io);

    for e in elem1'range loop
      put(io,"*SOLID SECTION, MATERIAL=MAT, ELSET=SET");
      put(io,e,0);
      new_line(io);
      put(io,elem1(e).area,2,6,3);
      new_line(io);
    end loop;

    put(io,"**");
    new_line(io);
    put(io,"** ANALYSIS STEP");
    new_line(io);
    put(io,"*STEP, AMPLITUDE=RAMP, NLGEOM=YES, INC=");
    put(io,neqit,0);
    new_line(io);

    put(io,"**");
    new_line(io);
    put(io,"** STATIC ANALYSIS");
    new_line(io);
    put(io,"* STATIC");
    new_line(io);
    put(io,(1.0 / long_float(nsubst)),2,6,3);
    put(io,",  1.000000E+00");
    new_line(io);

    put(io,"**");
    new_line(io);
    put(io,"** BOUNDARY CONDITIONS");
    new_line(io);
    put(io,"*BOUNDARY");
    new_line(io);

    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        if index(glob,axis) > n then
          put(io,glob,5);
          put(io,", ");
          put(io,axis,5);
          new_line(io);
        end if;
      end loop;
    end loop;

    put(io,"**");
    new_line(io);
    put(io,"** CONCENTRATED FORCES");
    new_line(io);
    put(io,"*CLOAD");
    new_line(io);

    for glob in 1 .. nn loop
      for axis in 1 .. 3 loop
        if force(glob,axis) /= 0.0 then
          put(io,glob,5);
          put(io,", ");
          put(io,axis,5);
          put(io,", ");
          put(io,force(glob,axis),2,6,3);
          new_line(io);
        end if;
      end loop;
    end loop;

    put(io,"**");
    new_line(io);
    put(io,"** PRINT DISPLACEMENT AND REACTION FORCE");
    new_line(io);
    put(io,"*NODE PRINT");
    new_line(io);
    put(io,"U");
    new_line(io);
    put(io,"RF");
    new_line(io);

    put(io,"**");
    new_line(io);
    put(io,"** PRINT STRESS");
    new_line(io);
    put(io,"*EL PRINT");
    new_line(io);
    put(io,"S");
    new_line(io);
    put(io,"*END STEP");
    new_line(io);

    close(io);
  end abaqus;

begin
  input;

  start;

  abaqus;

  qn.lbfgs(maxlist,n,maxminor,maxls,
           minor,
           maxnorm,lstol,
           f,
           x,g,
           optimal);

  autocad;

  output;
end iconvex2;
