///////////////////////////////////////////////////////////////
// Copyright G. Allaire, A. Karrman, G. Michailidis, January 2012
//
// A Scilab toolbox for 2-d structural optimization
// by the level set method with the topological gradient.
//
// Based on the work of G. Allaire, F. Jouve, A.-M. Toader,  
// J. Comp. Phys. Vol 194/1, pp.363-393 (2004). 
////////////////////////////////////////////////////////////////
//
// This file contains the parameters of the following test case:
// L-beam 
////////////////////////////////////////////////////////////////
//
clear ntop;
// PARAMETERS
nelx = 120 ;         // element count in x direction
nely = 120 ;         // element count in y direction
xlength = 1 ;       // working domain length
dx = xlength/nelx ; // x space step size
dy = dx ; // y space step size
yheight = dy*nely ;       // working domain height

hx = 7 ;            // x direction holes
hy = 5.5 ;            // y directino holes
r = .5 ;            // hole size (between 0 and 1)

eps = .001 ;        // "hole" material density
lagV = 150. ;        // the volume lagrangian multiplier
lagP = .0 ;        // the perimeter multiplier
e2 = 4*dx^2.;     // coefficient for the regularization in front of the Laplacian

nodex = linspace(0,xlength,nelx+1) ;  // x space for nodes
nodey = linspace(0, yheight,nely+1) ; // y space for nodes
FEx = linspace(0,xlength,nelx) ;      // x space for finite elements
FEy = linspace(0, yheight,nely) ;     // y space for finite elements

RIiterinit = 50 ; // reinitialization iterations when first reinitializing
RIiter = 5 ; // reinitialization iterations while solving HJ equation
RIfreq = 5 ;  // number of HJ iterations for each time we reinitialize
HJiter0 = 10 ; // original number of HJ iterations
entire = 50 ; // total number of optimization iterations
allow = .01 ; // fraction that the objective function can increase in shape advection
allow_top = .1 ; // fraction that the objective function can increase in topological derivation
ntop = 6 ; //number of advection steps between two topological gradient steps
percentage_in = 0.04; // percentage of the current volume to be removed in the first topological gradient step
percentage = 0.04; // percentage of the current volume to be removed in each topological gradient step

// FINITE ELEMENT MATRICES
KE = lk(); // the stiffness matrix
K = sparse([],[],[2*(nelx+1)*(nely+1), 2*(nelx+1)*(nely+1)]); // the global stiffness matrix
F = sparse([],[],[2*(nely+1)*(nelx+1),2]); // the matrix of applied forces
U = sparse([],[],[2*(nely+1)*(nelx+1),2]); // the vector displacement matrix

// FINITE DIFFERENCES MATRIX
K1 = sparse([],[],[(nelx+1)*(nely+1),(nelx+1)*(nely+1)]) ;  //We define the matrix of the velocity regularization

// SETTING OF THE ELASTICITY PROBLEM

  // THE FORCE
  // Each row corresponds to a different force.
  // The first two values are the fraction along the x and y axes of the working
  // domain (the origin is the top left corner).  The third value is binary indicating
  // whether the force is horizontal (0) or vertical (1).  The fourth value gives
  // the strength of the force and its direction (negative for leftward or downward
  // forces; positive for rightward or upward forces).
  forceMatrix = [.75 0 0 1] ;
  
  forceCount = size(forceMatrix,1) ; // We count the number of applied forces
  for force = 1 : forceCount
    F(c(forceMatrix(force,1:3)),force) = forceMatrix(force,4) ;
  end

  // FIXED BOUNDARIES WITH DIRICHLET CONDITIONS
  fixeddofs = [(nely+1):2*(nely+1)] ; // We fix x and y degrees of freedoms for nodes
  
  alldofs     = [1:2*(nely+1)*(nelx+1)];
  freedofs    = setdiff(alldofs,fixeddofs);
  
  // PASSIVE ELEMENTS
  // For this part we can make sure that certain areas in the working domain 
  // are either always part of the structure or always not part of the structure.
  function FEthetaOut = passive(FEthetaIn)
    FEthetaOut = FEthetaIn ;
    for i = 1 : ceil(nely*.5)
      for j = 1 : ceil(nelx*.5)
        FEthetaOut(i,j) = eps ;
      end
    end
  endfunction
  
//

// INITIALIZATION
phi0 = mesh0(hx, hy,r) ;