# Objective: convex quadratic
# Constraints: convex quadratic

########################################################################
# Continuum elements.
# Iterative procedure to impose upper bounds of 1 on sigmas.
########################################################################

    param lambda;
    param mu;
    param porosity;
    
    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 default {};
#    set ANCHORS within NODES := {};
#				{ x in X, y in Y : 
#			          x == 0 
##			          && y >= floor(m/3) && y <= m-floor(m/3) 
#			          };
    
    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 rho {ELEMENTS};
    param Vol, default 1;
    param reducedVol;

#################################################################
# nonlinear model for w, theta  (w is u and theta is xi in my notes)
#################################################################

    var w {NODES, 1..2};
    var theta >= 0;
    
    maximize work: 
        2 * sum {(x,y) in NODES, d in 1..2} load[x,y,d]*w[x,y,d]
	- sum {(x,y) in ELEMENTS: rho[x,y] == 1}
		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] 
		)
	- theta
	;
    
    subject to element_eqs {(x,y) in ELEMENTS: rho[x,y] != 1}: 
        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] 
        )
        <= theta/reducedVol;
    
    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}:
#    subject to compat {(x,y) in NODES, d in 1..2}: 
         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.95;
    
#    let n := 36;
#    let m := 24;	# must be even
    let n := 18;
    let m := 18;	# must be even

#    let load[n,floor(m/2),2] := -1;
    let {(x,y) in NODES: x==0 && y>0 && y<m} load[x,y,1] := -y*(m-y);
    let {(x,y) in NODES: x==n && y>0 && y<m} load[x,y,1] :=  y*(m-y);
    let {(x,y) in NODES: y==0 && x>0 && x<n} load[x,y,2] := -x*(n-x);
    let {(x,y) in NODES: y==m && x>0 && x<n} load[x,y,2] :=  x*(n-x);
#    let load[0,0,1] := -1/2; let load[0,0,2] := -1/2;
#    let load[n,0,1] :=  1/2; let load[n,0,2] := -1/2;
#    let load[n,m,1] :=  1/2; let load[n,m,2] :=  1/2;
#    let load[0,m,1] := -1/2; let load[0,m,2] :=  1/2;

    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, theta, work, element_eqs; #, anchor;
    option solver loqo;
    option loqo_options "verbose=2 inftol=1.0e-3 honor_initsol pred_corr=0";
    option presolve 1;

    problem linmodel: sigma, zero, compat, tot_vol;
    option solver loqo;
    option loqo_options "verbose=2 sigfig=0 inftol=1.0e-5 honor_initsol convex";
    option presolve 0;

    let { (x,y) in NODES, d in 1..2 } w[x,y,d] := 0.001;
    let { (x,y) in ELEMENTS } rho[x,y] := 0.5;
    let reducedVol := Vol;
    
    param iter;
    let iter := 1;
    repeat {
	printf: "ITERATION: %d \n", iter;
        solve nlmodel;

        let { (x,y) in NODES, d in 1..2 } u[x,y,d] := w[x,y,d];
        let { (x,y) in ELEMENTS: rho[x,y] == 1 } sigma[x,y] := 1;
	fix { (x,y) in ELEMENTS: rho[x,y] == 1 } sigma[x,y];
        
        solve linmodel;
        
        let {(x,y) in ELEMENTS} 
	    rho[x,y] := if (sigma[x,y] > 0.9999) then 1 else sigma[x,y];
	let reducedVol := Vol - sum {(x,y) in ELEMENTS: rho[x,y] == 1} 1 ;

	display max {(x,y) in ELEMENTS} sigma[x,y];

	let iter := iter+1;
    } while (max {(x,y) in ELEMENTS} sigma[x,y] > 1.01);
    
    option display_eps 0.0001;
    printf: "%d \n", 
    card({(x,y) in ELEMENTS: rho[x,y] > 1.0e-8})
    > structure5.out;
    printf {(x,y) in ELEMENTS: rho[x,y] > 1.0e-8}: 
    "%3d %3d %10.4f \n", x, y, rho[x,y]
    > structure5.out;