# Objective: convex quadratic
# Constraints: convex quadratic

########################################################################
# Here we just compute the compliance for some specific settings for sigma.
#
# Continuum elements.
########################################################################

    param lambda;
    param mu;
    param porosity;
    param length;
    
    param kk;	# each super-element is a kxk block of elements
    param n;	# number of elements across x direction
    param m;	# number of elements across y direction
    
    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 {};
    
    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 && rho[x,y] !=0 }: 
        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;
    
#################################################################
# linear model for compliance
#################################################################

    minimize compliance: 
	sum{(x,y) in NODES, d in 1..2} load[x,y,d]*w[x,y,d]*length^2;
    
    subject to compat2 {(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]*w[x+k-i,y+l-j,1] +
    	 K[i,j,k,l,d,2]*w[x+k-i,y+l-j,2] 
         )*rho[x-i,y-j]
         = load[x,y,d];
    
#################################################################
# set data
#################################################################

    let lambda := 1;
    let mu := 1;
    let porosity := 0.50;
    let length := 1;
    
    let kk := 2;
    let n  := 32;
    let m  := n;	

    let {(x,y) in NODES: x==0 && y mod kk==0} load[x,y,1] := -kk;
    let {(x,y) in NODES: x==n && y mod kk==0} load[x,y,1] :=  kk;
    let {(x,y) in NODES: y==0 && x mod kk==0} load[x,y,2] := -kk;
    let {(x,y) in NODES: y==m && x mod kk==0} load[x,y,2] :=  kk;
    let load[0,m,1] := 0; let load[0,m,2] := 0;
    let load[n,0,1] := 0; let load[n,0,2] := 0;

    display (32/n) * sum {(x,y) in NODES, d in 1..2} abs(load[x,y,d]);

    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 sigfig=7";
    option presolve 1;

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

    problem compliance_mod: w, compliance, compat2;
    option solver loqo;
    option loqo_options "verbose=2 inftol=1.0e-3 sigfig=7";
    option presolve 1;

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

#        let { (x,y) in ELEMENTS: (floor(x/kk)+floor(y/kk)) mod 2 == 0 } 
#		rho[x,y] := 1.0;

        let { (x,y) in ELEMENTS: 
	    (floor(x/kk)-(floor(n/kk)-1)/2)^2 +
	    (floor(y/kk)-(floor(m/kk)-1)/2)^2 >= 43 } rho[x,y] := 1.0;

#        let { (x,y) in ELEMENTS } rho[x,y] := 8/9;

        display card({(x,y) in ELEMENTS: rho[x,y] == 0});

        option display_eps 0.0001;
        printf: "%d \n", 
        card({(x,y) in ELEMENTS: rho[x,y] > 1.0e-8})
        > structure9.out;
        printf {(x,y) in ELEMENTS: rho[x,y] > 1.0e-8}: 
        "%3d %3d %10.4f \n", x, y, rho[x,y]
        > structure9.out;
    
        solve compliance_mod;