# Torsion problem
# Liz Dolan - Summer 2000
# Version 2.0 - October 2000

model;

param nx > 0, integer; # grid points in 1st direction
param ny > 0, integer; # grid points in 2nd direction
param c;               # constant

param hx := 1/(nx+1);     # grid spacing
param hy := 1/(ny+1);     # grid spacing
param area := 0.5*hx*hy;  # area of triangle

# Distance to the boundary.

param D {i in 0..nx+1,j in 0..ny+1} := 
  min(min(i,nx-i+1)*hx,min(j,ny-j+1)*hy);

# v defines the finite element approximation.

var v {i in 0..nx+1,j in 0..ny+1};

var linLower = sum {i in 0..nx,j in 0..ny} (v[i+1,j] + v[i,j] + v[i,j+1]);

var linUpper = sum {i in 1..nx+1,j in 1..ny+1} (v[i,j] + v[i-1,j] + v[i,j-1]);

var quadLower = sum {i in 0..nx,j in 0..ny} (
  ((v[i+1,j] - v[i,j])/hx)**2 + ((v[i,j+1] - v[i,j])/hy)**2);

var quadUpper = sum {i in 1..nx+1, j in 1..ny+1}(
  ((v[i,j] - v[i-1,j])/hx)**2 + ((v[i,j] - v[i,j-1])/hy)**2);

minimize Stress:
  area*((quadLower + quadUpper)/2 - c*((linLower + linUpper)/3.0));

subject to distanceBound {i in 0..nx+1,j in 0..ny+1}:
  -D[i,j] <= v[i,j] <= D[i,j];

# Torsion problem
# Liz Dolan
# Version 2.0 - October 2000

data;

# Set the number of discretization points and design parameters

param nx := 600;
param ny := 600;
param c := 5;

# Set the starting point to the upper bound.

let {i in 0..nx+1,j in 0..ny+1} v[i,j]:= D[i,j];