%Clay_SSL
clear
clf
L=10;
N=100; 
xx = 0:(L/(N-1)):L;
yy = 0:(L/(N-1)):L;
f = 0 ;

%Initial Conditions
S=ones(N,N);            %No Suction anywhere
q=.5;                   %Evaporation rate (Suction flux)
dy = L/(N-1);
Scrit = 35;           %Critical stress to form a crack
S = rand(N,N)*10;       %add noise

S(1:N,1)= 0;           %No suction on bottom 
%  S(1,1:N) = 0;       %No Suction on right side
%  S(N,1:N)= 0;        %No Suction on left side



for iter= (1:515);
    Snew = S;
    for j=(2:N-1);  
        for k=(2:N-1);
            if S(j,k) ~= 0  
                Snew(j,k)=0.25*(S(j-1,k) - 4*S(j,k) + S(j+1,k)+ S(j,k+1)+S(j,k-1)) + S(j,k); 
            end
        end
    end
%     Make interior of crack a flux area 
    for j = (3:N-2)
        for k = (3:N-2)
            if S(j+2,k)==0||S(j-2,k)==0||S(j,k+2) == 0 && S(j,k) ~= 0
               Snew(j,k) = 0.25*(S(j-1,k) - 4*S(j,k) + S(j+1,k)+ S(j,k+1)+S(j,k-1)) + S(j,k);
            end
        end
    end
    S=Snew;
%     Form a crack once it reaches a critical stress level Scrit
     for j = (2:N-1);
         for k = (2:N); 
             if S(j,k) > Scrit 
                S(j,k)  = 0;   
                S(j-1,k) = .5*S(j-1,k);
                S(j+1,k) = .5*S(j+1,k);
                S(j,k-1)=0;
             end             
             if S(j,k) == 0
                S(j,k-1) = 0;
             end
         end
             
     end
         for k = (2:N);
             if S(j+1,k) == 0 || S(j-1,k) == 0  && S(j,k)~=0;       %|If any of my neighbors are 0,become flux
                S(j,k) =  S(j,k) + q*dy; 
             end
      
         end
      

     
%Flux Boundary Conditions at the top (evaporation)
     for j = (1:N);
         if S(j,N) ~= 0
            S(j,N) = q*dy + S(j,N);
         end
     end 
         
       %No flux boundary conditions         
        S(1,1:N) = S(2,1:N);            %Left Side        
        S(N,1:N)= S(N-1,1:N);           %Right Side

       %Plot commands
         [X,Y] = meshgrid(yy,xx);
         if mod(iter+9,10) == 0
             f = f+1;
             surf(X,Y,S)
              saveas(gcf,['DCNSO',sprintf('%03d',f),'.png'])
             pause (.1)
         end
         view([90 -90]);
end
 

 S = S';
