$onText
Determine the stress potential in an infinitely long cylinder when torsion
is applied.

This model is from the COPS benchmarking suite.
See http://www-unix.mcs.anl.gov/~more/cops/.

The number of internal grid points can be specified using the command
line parameters --nx and --ny.
COPS performance tests have been reported for nx-1 = 50,
ny-1 = 25, 50, 75, 100

References:
Dolan, E D, and More, J J, Benchmarking Optimization Software with COPS.
  Tech. rep., Mathematics and Computer Science Division, 2000.
Glowinski, R, Numerical Methods for Nonlinear Variational Problems.
  Springer Verlag, 1984.
Averick, B M, Carter, R G, More, J J, and Xue, G L, The MINPACK-2 Test
  Problem Collection. Tech. rep., Mathematics and Computer Science
  Division, Argonne National Laboratory, 1992.
Andrei, N., Criticism of the unconstrained optimization algorithms reasoning.
  Romanian Academic Press, Bucharest, 2009, page: 432-434.
$offText

$if not set nx $set nx 100
$if not set ny $set ny 100

sets nx grid points in 1st direction / x0*x%nx% /
     ny grid points in 2st direction / y0*y%ny% /

alias(nx,i),(ny,j);

parameters D(nx,ny) Distance to the boundary
           hx grid spacing for x
           hy grid spacing for y
           area area of triangle
           c some constant / 5.0 /;

hx := 1/(card(nx)-1);
hy := 1/(card(ny)-1);
area := 0.5*hx*hy;

D(i,j) := min(min(ord(i)-1,card(nx)-ord(i))*hx,
              min(ord(j)-1,card(ny)-ord(j))*hy);

variables v(nx,ny) the finite element approximation stress,
          linLower,
          linUpper,
          quadLower,
          quadUpper,
          stress;

Equations defLL,
          defLU,
          defQL,
          defQU,
          defstress;

defLL.. linLower =e= sum((nx(i+1),ny(j+1)), v[i+1,j] + v[i,j] +
                                            v[i,j+1]);

defLU.. linUpper =e= sum((nx(i-1),ny(j-1)), v[i,j] + v[i-1,j] +
                                            v[i,j-1]);

defQL.. quadLower =e= sum((nx(i+1),ny(j+1)),
                                    sqr((v[i+1,j]-v[i,j])/hx) +
                                    sqr((v[i,j+1]-v[i,j])/hy));

defQU.. quadUpper =e= sum((nx(i-1),ny(j-1)),
                                    sqr((v[i,j]-v[i-1,j])/hx) +
                                    sqr((v[i,j]-v[i,j-1])/hy));

defstress.. stress =e= area*( (quadLower + quadUpper)/2 -
                          c*(linLower + linUpper )/3);

model torsion / all /;

$ifThenI x%mode%==xbook
v.lo(i,j) = -d(i,j);
v.up(i,j) = d(i,j);
v.l (i,j) = d(i,j);

display d,hx,hy,area;
$onEcho >minos.opt
  superbasics limit = 5000
$offEcho
torsion.workspace=120;
$endIf

solve torsion minimizing stress using nlp;

$ifThenI x%mode%==xbook
file res1 /torsion.txt/;
res1.pw=4000;
put res1;

loop(nx, loop(ny, put v.l(nx,ny):6:2; ); put /;); put /;
$endIf

* End Torsion