////////////////////////////////////////////////////////////////
// Copyright G. Allaire, A. Karrman, G. Michailidis, January 2012
//
// 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 is the main routine to execute in Scilab.
// Before that you should have load the subroutines 
// in the file functions.sci as well as the parameters 
// of the desired test case.
// This program performs compliance minimization with 
// constraints on the volume and perimeter. These constraints 
// are handled by fixed Lagrange multipliers.
// 
////////////////////////////////////////////////////////////////
//
// Check if the variable 'ntop' for topological gradient is defined.
if isdef('ntop') then
  ntop = ntop;
else
  ntop = 0;
end
// REINITIALIZATION
phi00 = mesh00(phi0,RIiterinit) ;

// Set phi to the reinitialized state:
phi = phi00 ;

// Plot the initialization:
scf(0);
clf()
xset('colormap',graycolormap(10))
title('Initialization')
grayplot(FEx,-FEy,-passive(FEdensity(phi))',axesflag = 2) ;

filename='Initialization';
xs2jpg(0,filename);

printf('\nOptimization started\n') ;
stacksize('max') ; // We want to make sure that we have enough memory available (only with Scilab version 5.1.1).

// FE ANALYSIS

// Define the elements' densities based on phi.
FEtheta = FEdensity(phi,eps) ;   
// Set unchanging densities depending on the 'passive' function.
// This 'passive' zone is defined in the parameter file of the 
// considered test case. 
FEtheta = passive(FEtheta) ;  

if ntop==0 then
  // shape derivative
  // The output of the finite element analysis is the elastic energy density:
  // lvlAeueu is that field defined on nodes (it is used as the velocity in
  // transport level set equation), FEAeueu is the same field defined on 
  // elements (it is used to compute the compliance).
  [lvlAeueu,FEAeueu] = FE(FEtheta,KE,K,F,U) ; 
else
  [FEAeueu,FEAeueucomp] = FEtop(FEtheta,KE,K,F,U);
  // topological derivative
  // In this case, FEAeueu denotes the topological gradient,
  // while FEAeueucomp contains the energy density defined on elements.
end
// Define the velocity field:
if ntop==0 then
  // shape gradient
  V = lvlAeueu/(dx*dy) - lagV;
  // regularization of the velocity field
  V = regularize(phi,V,e2,K1,lagP);
else
  // topological gradient
  V = FEAeueu - lagV;
end


// CALCULATE THE OBJECTIVE FUNCTION
if ntop==0 then
  totcomp = compliance(FEAeueu);
else
  totcomp = compliance(FEAeueucomp) ;
end
totvol = volume(FEtheta) ;
objective = lagV*totvol+lagP*perimeter(phi)+totcomp ;

// We track the objective function after each optimization iteration using 'objectivePlot'
objectivePlot = objective ;


// OPTIMIZATION LOOP

i = 1 ; // The current optimization iteration 
HJa = 1 ; // The level set solution "attempt"
HJiter = HJiter0 ; // The initial number of time steps for solving the transport level set equation
e3 = 1;  // Coefficient used to reduce, if neccessary, the advection time step below the 
// limit imposed by the cfl condition
allow_adv = allow;

while i<=entire  
  
  if (ntop~=0) then
    if modulo(i,ntop)==0 | (i==1) then
      // Test shape using the topological derivative:
      if (i==1) then 
        phiTest = solvelvlset_top(phi,V,percentage_in,FEtheta);
      else     
        phiTest = solvelvlset_top(phi,V,percentage,FEtheta);
      end
    else
      dt = 0.5*e3*min(dx,dy)/max(abs(V)) ;
      phiTest = solvelvlset(phi,V,dt,HJiter,lagP) ;  
    end
  else
    // Solve the transport level set equation using the velocity V:
    dt = 0.5*e3*min(dx,dy)/max(abs(V)) ;
    phiTest = solvelvlset(phi,V,dt,HJiter,lagP) ;
  end
  
  // FE ANALYSIS 
  // Define material density based on phi
  FEthetaTest = FEdensity(phiTest,eps) ; 
  // Set passive element densities to 'eps'
  FEthetaTest = passive(FEthetaTest) ;
  // Perform the finite element analysis and compute 
  // the elastic energy density:
  if (ntop~=0) then
    if modulo(i,ntop)==0 | (i==1) then
      [FEAeueuTest,FEAeueucompTest] = FEtop(FEthetaTest,KE,K,F,U);
    else
      [lvlAeueuTest,FEAeueuTest] = FE(FEthetaTest,KE,K,F,U); 
    end
  else
    [lvlAeueuTest,FEAeueuTest] = FE(FEthetaTest,KE,K,F,U);
  end

  // CALCULATE THE OBJECTIVE FUNCTION
  if (ntop~=0) then
    if modulo(i,ntop)==0 | (i==1) then
      totcompTest = compliance(FEAeueucompTest) ;
    else
      totcompTest = compliance(FEAeueuTest) ;
    end
  else
      totcompTest = compliance(FEAeueuTest) ;
  end
  totvolTest = volume(FEthetaTest) ;
  objectiveTest = lagV*totvolTest+lagP*perimeter(phiTest)+totcompTest ;
  // Plot the test shape: 
  scf(0) 
  clf()
  xset('colormap',graycolormap(10)) ;
  title('Test shape')
  grayplot(FEx,-FEy,-FEthetaTest',axesflag = 2) ;
  printf('\nIteration %d of %d,  HJ attempt %d, HJiter %d, objective = %f, objectiveTest = %f, volume = %f\n',...
    i,entire,HJa,HJiter,objective,objectiveTest,totvolTest) ;
    
  // OBJECTIVE FUNCTION MUST DECREASE (up to some tolerance 'allow')
  if (ntop~=0) then
    if modulo(i,ntop)==0 | (i==1) then
      allow = allow_top;
    else
      allow = allow_adv;
    end
  else
    allow = allow_adv;
  end 
  if i>=(entire*3/4) then
    allow = 0;  // After some iterations, we switch-off the 'allow' parameter,
                // in order to converge.
    ntop =0;
  end
  if objectiveTest <= objective*(1+allow) then
    // The current design is OK: move on to the next iteration, 
    // using the test versions as our new variables to work with
    i=i+1; // Move on to the next optimization iteration
    phi = phiTest ;
    FEtheta = FEthetaTest ;
    objective = objectiveTest ;
    if (ntop~=0) then
      if modulo(i,ntop)==0 | modulo(i-1,ntop)==0 | i==2 then
        if modulo(i,ntop)==0 then
          [FEAeueu,FEAeueucomp] = FEtop(FEtheta,KE,K,F,U);
          V = FEAeueu - lagV;
        else
          [lvlAeueu,FEAeueu] = FE(FEtheta,KE,K,F,U);
          V = lvlAeueu/(dx*dy) - lagV;
          V = regularize(phi,V,e2,K1,lagP);
        end
      else
        lvlAeueu = lvlAeueuTest ;
        FEAeueu =FEAeueuTest ;
        V = lvlAeueu/(dx*dy) - lagV;
        V = regularize(phi,V,e2,K1,lagP);
      end
    else
       lvlAeueu = lvlAeueuTest ;
       FEAeueu =FEAeueuTest ;
       V = lvlAeueu/(dx*dy) - lagV;
       V = regularize(phi,V,e2,K1,lagP);
    end
    HJiter = min(10,max(floor(HJiter*1.1),HJiter+1)) ;
    objectivePlot($+1) = objective ;
    // Plot the new shape:  
    clf()
    xset('colormap',graycolormap(10)) ;
    title('Evolving shape')
    grayplot(FEx,-FEy,-FEtheta',axesflag = 2) ;
    /////////////////////////////////////
    HJa = 1 ; // Reset the lvl-set "attempt" variable
    e3 = 1;
  else
    i=i+1;
    if (ntop~=0) then
      if modulo(i,ntop)==0 | modulo(i-1,ntop)==0 then
        if modulo(i,ntop)==0 then
          [FEAeueu,FEAeueucomp] = FEtop(FEtheta,KE,K,F,U);
          V = FEAeueu - lagV;
        else 
          [lvlAeueu,FEAeueu] = FE(FEtheta,KE,K,F,U);
          V = lvlAeueu/(dx*dy) - lagV;
          V = regularize(phi,V,e2,K1,lagP);
        end
      else
        // The current design is bad: try again, this time with fewer 
        // time steps for the transport level set equation
        HJiter = floor(HJiter/2) ;
        if HJiter == 0 then
           HJiter = 1 ;
           e3 = e3/2;
        end    
        HJa = HJa + 1 ; // Increment the lvl-set "attempt" variable
      end
    else
      // The current design is bad: try again, this time with fewer 
      // time steps for the transport level set equation
      HJiter = floor(HJiter/2) ;
      if HJiter == 0 then
         HJiter = 1 ;
         e3 = e3/2;
      end    
      HJa = HJa + 1 ; // Increment the lvl-set "attempt" variable
    end
  end
end

// Plot the final shape:  
clf()
xset('colormap',graycolormap(10)) ;
title('Final shape')
grayplot(FEx,-FEy,-FEtheta',axesflag = 2) ;
scf(0);
filename='Finaldesign';
xs2jpg(0,filename);
printf('\n Export of Final Design\n') ;

// Plot the objective function:
clf()

xtitle('Convergence history','Iteration','Objective function')
plot((1:length(objectivePlot)) - 1,objectivePlot) ;

printf('\nOptimization finished\n') ;