///////////////////////////////////////////////////////////////
// Copyright G. Allaire, A. Karrman, October 2009
//
// A Scilab toolbox for 2-d structural optimization
// by the level set method.
//
// 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:
// bridge 
////////////////////////////////////////////////////////////////
//
clear ntop;
// PARAMETERS
nelx = 80 ;         // element count in x direction
nely = 48 ;         // element count in y direction
xlength = 2 ;       // working domain length
dx = xlength/nelx ; // x space step size
dy = dx ; // y space step size
yheight = dy*nely ;       // working domain height

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

eps = .001 ;        // "hole" material density
lagV = 15. ;       // the volume Lagrange multiplier
lagP = 0.01 ;          // the perimeter Lagrange 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 ; // number of time steps in the reinitialization of the initial design
RIiter = 5 ; // number of time steps in the reinitialization of further designs
RIfreq = 5 ;  // frequency of reinitialization, i.e. number of time steps in the 
// transport level set equation between two reinitializations
HJiter0 = 30 ; // original number of transport time steps 
entire = 30 ; // total number of optimization iterations
allow = .02 ; // fraction that the objective function can increase

// 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 = [.5 1 1 -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 = [c([0 1 1]) c([0 1 0]) c([1 1 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 ;
  endfunction
  
//

// INITIALIZATION
phi0 = mesh0(hx, hy,r) ;
phi0(($-floor(.08*nely)):$,:) = -.1 ;