# Objective: linear
# Constraints: convex quadratic

########################################################################
# Continuum elements.
# Works fine with presolve turned on.  Fails for some reason when 
# presolve turned off.
########################################################################

param lambda;
param mu;
param porosity;

param n;
param m;	# must be even
#param n;
#param m;	# must be even

param N := if (n>m) then n else m;

set X0 := {0..n-1};
set Y0 := {0..m-1};

set X := {0..n};
set Y := {0..m};

set NODES := X cross Y;		# A lattice of Nodes

set ANCHORS within NODES := { x in X, y in Y : 
			      x == 0 
			      && y >= floor(m/3) && y <= m-floor(m/3) 
			      };
#display ANCHORS;

param load {(x,y) in NODES, d in 1..2} default 0;

set ELEMENTS := X0 cross Y0;

param xi {i in 0..1} := if (i==0) then -1 else 1;

param axx {i in 0..1, j in 0..1, k in 0..1, l in 0..1} 
    := xi[i]*xi[k]*(1 + xi[j]*xi[l]/3);
param ayy {i in 0..1, j in 0..1, k in 0..1, l in 0..1} 
    := xi[j]*xi[l]*(1 + xi[i]*xi[k]/3);
param axy {i in 0..1, j in 0..1, k in 0..1, l in 0..1} 
    := xi[i]*xi[l];
param ayx {i in 0..1, j in 0..1, k in 0..1, l in 0..1} 
    := xi[j]*xi[k];

param K {0..1, 0..1, 0..1, 0..1, 1..2, 1..2};
param u {NODES, 1..2};
param Vol, default 1;

#################################################################
# nonlinear model for w (w is u in my notes)
#################################################################

var w {NODES, 1..2};

maximize work: 
    sum {(x,y) in NODES, d in 1..2} load[x,y,d]*w[x,y,d];

subject to element_eqs {(x,y) in ELEMENTS}: 
    sum {i in 0..1, j in 0..1, k in 0..1, l in 0..1}
    (
	w[x+i,y+j,1]*K[i,j,k,l,1,1]*w[x+k,y+l,1] +
	w[x+i,y+j,2]*K[i,j,k,l,2,1]*w[x+k,y+l,1] +
	w[x+i,y+j,1]*K[i,j,k,l,1,2]*w[x+k,y+l,2] +
	w[x+i,y+j,2]*K[i,j,k,l,2,2]*w[x+k,y+l,2] 
    )
    <= 1/Vol;

subject to anchor {(x,y) in ANCHORS, d in 1..2}:
    w[x,y,d] = 0;

#################################################################
# linear model for sigma
#################################################################

var sigma {ELEMENTS} >= 0;

maximize zero: 0;

subject to compat {(x,y) in NODES, d in 1..2: (x,y) not in ANCHORS}:
     sum {i in 0..1, j in 0..1, k in 0..1, l in 0..1: (x-i,y-j) in ELEMENTS}
     (
	 K[i,j,k,l,d,1]*u[x+k-i,y+l-j,1] +
	 K[i,j,k,l,d,2]*u[x+k-i,y+l-j,2] 
     )*sigma[x-i,y-j]
     = load[x,y,d];

subject to tot_vol:
    sum { (x,y) in ELEMENTS } sigma[x,y] = Vol;

#################################################################
# set data
#################################################################

let lambda := 1;
let mu := 1;
let porosity := 0.5;

#let n := 9;
#let m := 5;	# must be even
#let n := 18;
#let m := 10;	# must be even
let n := 36;
let m := 20;	# must be even
#let n := 72;
#let m := 40;	# must be even

let load[n,floor(m/2),2] := -1;

let Vol := porosity*n*m;

let {i in 0..1, j in 0..1, k in 0..1, l in 0..1}
    K[i,j,k,l,1,1] := (lambda + 2*mu)*axx[i,j,k,l] + mu*ayy[i,j,k,l];
let {i in 0..1, j in 0..1, k in 0..1, l in 0..1}
    K[i,j,k,l,2,2] := (lambda + 2*mu)*ayy[i,j,k,l] + mu*axx[i,j,k,l];
let {i in 0..1, j in 0..1, k in 0..1, l in 0..1}
    K[i,j,k,l,1,2] := lambda*axy[i,j,k,l] + mu*ayx[i,j,k,l];
let {i in 0..1, j in 0..1, k in 0..1, l in 0..1}
    K[i,j,k,l,2,1] := lambda*ayx[i,j,k,l] + mu*axy[i,j,k,l];

#################################################################
# control logic
#################################################################

problem nlmodel: w, work, element_eqs, anchor;
#option loqo_options "timing=1 verbose=2 inftol=1.0e-3 \
#	honor_initsol pred_corr=0";
#option minos_options "timing=1 superbasics=2000";
#option lancelot_options "timing=1 maxit=10000 outlev=0";
#option solver loqo;
#option solver minos;
#option solver lancelot;
option presolve 1;

problem linmodel: sigma, zero, compat, tot_vol;
#option solver loqo;
#option solver minos;
#option solver lancelot;
option loqo_options "verbose=2 sigfig=0"; 
option presolve 0;

let {(x,y) in NODES, d in 1..2} w[x,y,d] := 0.001;

solve nlmodel;

let {(x,y) in NODES, d in 1..2} u[x,y,d] := work * w[x,y,d];
#let {(x,y) in NODES, d in 1..2} u[x,y,d] := w[x,y,d];

solve linmodel;

option display_eps 0.0;

printf: "%d \n", 
    card({(x,y) in ELEMENTS: sigma[x,y] > 1.0e-8})
    > structure4.out;
printf {(x,y) in ELEMENTS: sigma[x,y] > 1.0e-8}: 
    "%3d %3d %10.4f \n", x, y, sigma[x,y]
    > structure4.out;

printf {x in X0, y in Y0}: "%10.6f \n",
    sigma[x,y] 
    > thickness.wrl;
printf {x in X0, y in Y0}: "%10.6f \n",
    -sigma[x,y] 
    > mthickness.wrl;

#printf {x in X, y in Y}: "%10.6f \n",
#    (sum {i in 0..1, j in 0..1: x-i in X0 && y-j in Y0} sigma[x-i,y-j])
#    /
#    card{i in 0..1, j in 0..1: x-i in X0 && y-j in Y0}
#    > thickness.wrl;
#printf {x in X, y in Y}: "%10.6f \n",
#    -(sum {i in 0..1, j in 0..1: x-i in X0 && y-j in Y0} sigma[x-i,y-j])
#    /
#    card{i in 0..1, j in 0..1: x-i in X0 && y-j in Y0}
#    > mthickness.wrl;
#display sigma;