function [Ux_All, Uy_All, M_All, F_All, Time_All, PosX, PosY, Front_All, FrontVelocity_All, Global_Area_All,Area_All,InitialMass,str] = GetSolutionStokesActin(config,SelectedTime)

global V0 M0 K A0 FX FY Actin_Velocity NormalX NormalY TotalArea  DomainX DomainY AreaExponent;

ticID = tic;

M_tot                = config.M_tot;
V0                   = config.V0;
M0                   = config.M0;
K                    = config.K;
A0                   = config.A0;
Domain               = config.Domain;
D                    = config.D;                        
eta                  = config.eta;                        
ksi                  = config.ksi;                          
sigma                = config.sigma;
zeta                 = config.zeta;
epsilon              = config.epsilon;
koff                 = config.koff;
F0                   = config.F0;

InitCondition        = config.InitCondition;
T_max                = config.T_max; 
MeshSize             = config.MeshSize;
Front.Formula        = config.FrontFormula;
NumFrames            = config.NumFrames;
ExactPolygon         = [];

MaxMyosinAdvection   = config.MaxMyosinAdvection;
MaxMyosinDiffusion   = config.MaxMyosinDiffusion;
AreaExponent         = config.AreaExponent;

M_tot_input = M_tot;

DomainX = Domain{1};
DomainY = Domain{2};
NumNodesX = round((DomainX(2)-DomainX(1))/MeshSize);
NumNodesY = round((DomainY(2)-DomainY(1))/MeshSize);
NumNodes  = NumNodesX*NumNodesY;
PosX = repmat((1:NumNodesX)', 1, NumNodesY);
PosX = DomainX(1)+(PosX(:) - 1/2)*MeshSize;
PosY = repmat(1:NumNodesY, NumNodesX,1);
PosY = DomainY(1)+(PosY(:) - 1/2)*MeshSize;

% Set time step size  
MaxVelocity = 1;
dt = config.TimeStepScale*MeshSize/(MaxVelocity*10);
NumTimeNodes = round(T_max/dt);
dt = T_max/NumTimeNodes; 
NumTimeNodes=round(T_max/dt);

% modify for precise last time output
det_last_output = 0;
iT_last_output = 0;
if  mod(NumTimeNodes,max(round(NumTimeNodes/NumFrames),1))>0
  iT_last_output = NumTimeNodes;
  det_last_output = 1;
  NumFramesFinal = floor(NumTimeNodes/max(round(NumTimeNodes/NumFrames),1))+2;
else
  NumFramesFinal = NumTimeNodes/max(round(NumTimeNodes/NumFrames),1)+1;
end

% file name
if config.FixedBoundary
    str = [config.VelBoundaryCond '_' config.ActinBoundaryCond '_fxd_N=' num2str(config.NumNodesX) '_T=' ...
        num2str(T_max) '_V0=' num2str(V0) '_M_tot=' num2str(M_tot_input) ...
        '_eta=' num2str(eta) '_eps=' num2str(epsilon) '_zeta=' num2str(zeta) ...
        '_D=' num2str(D) '_koff=' num2str(koff) '_F0=' num2str(F0) '_' config.FrontRedistMethod];   
else
    str = [config.VelBoundaryCond '_' config.ActinBoundaryCond '_N=' num2str(config.NumNodesX) '_T=' ...
        num2str(T_max) '_V0=' num2str(V0) '_M_tot=' num2str(M_tot_input) ...
        '_eta=' num2str(eta) '_eps=' num2str(epsilon) '_zeta=' num2str(zeta) ...
        '_D=' num2str(D) '_koff=' num2str(koff) '_F0=' num2str(F0) '_' config.FrontRedistMethod];
    
end
   
if config.AvoidSingbyDiffu 
 str = [str '_' config.DiffCoefApproxMeth '-' num2str(MaxMyosinDiffusion)];
end

if config.AvoidSingbyAdvec
 str = [str '_ad-' num2str(MaxMyosinAdvection)];
end

CurrentDir = pwd;
if ispc
    delim_ch = '\'; 
else
    delim_ch = '/'; 
end
if ~exist('actin', 'dir')
    mkdir([CurrentDir delim_ch 'actin' delim_ch]);
end

% preparing front velocity at t=0
% different formulas for vp in ZS model
if config.bilinearInterp
    Actin_Velocity = @(x,y) [InitCondition.Ux(x,y) InitCondition.Uy(x,y)];          
else
    FX = @(x,y) InitCondition.Ux(x,y);
    FY = @(x,y) InitCondition.Uy(x,y);
end

switch config.InitDomainShape
  case 'Circle'
    TotalArea = config.InitArea;
    NormalX   = @(x,y) x/norm([x y]);
    NormalY   = @(x,y) y/norm([x y]);
  case 'ShapeFromSavedData'
    TotalArea = config.InitArea;
    iNeigh = 1;
    InitFront = config.InitFront;
    Tangents  = [InitFront.Polygon(end-iNeigh+1:end,:); InitFront.Polygon(1:end-iNeigh,:)] - [InitFront.Polygon(1+iNeigh:end,:); InitFront.Polygon(1:iNeigh,:)];
    Normals   = [-Tangents(:,2), Tangents(:,1)]./repmat(sqrt(Tangents(:,1).^2+Tangents(:,2).^2),1,2);
    NormalX   = scatteredInterpolant(InitFront.Polygon(:,1), InitFront.Polygon(:,2), Normals(:,1),'nearest');
    NormalY   = scatteredInterpolant(InitFront.Polygon(:,1), InitFront.Polygon(:,2), Normals(:,2),'nearest'); 
end


% Initialize front, mesh & solution     

[f, Front] = InitializeFront(Front, ExactPolygon, Domain, ...
  round(config.FrontRefinementScale*NumNodesX), T_max, config.FrontPropagate, ...
  config.FrontRedistMethod, config.FrontRedistFrequency);
if isempty(Front.Polygon)
    disp('Front variable is empty. Please choose a smaller MeshSize.')
    Ux_All=[]; Uy_All=[]; M_All=[]; F_All=[]; Time_All=[]; PosX=[]; PosY=[]; Front_All=[]; FrontVelocity_All=[]; Global_Area_All=[];
    return;
end

[Ux, Uy, M, F, InsideOldMask, InteriorOldMask, BoundaryOldMask] =...
    InitializeSolutionActin(PosX, PosY, Front, InitCondition, NumNodesX);
[NatNeighborPairs, NatNeighborBoundaryPairs, Vertices, VoronoiCells, Global_Area, D_Inv, Directions, MidPoints] =...
    InitializeControlVolumes(PosX, PosY, InsideOldMask, InteriorOldMask, BoundaryOldMask, Front, MeshSize, NumNodesX, NumNodesY);

InitialMass = sum(M.*Global_Area);
% making sure that the initial solution has the correct M_total
if strcmp(config.InitialCondType,'SolutionFromSavedData')
  M = config.M_tot/InitialMass*M;
end
  
% modify for precise last time output
if config.returnAllData
    Ux_All = sparse(NumNodes, NumFramesFinal);
    Ux_All(:,1) = Ux;
    Uy_All = sparse(NumNodes, NumFramesFinal);
    Uy_All(:,1) = Uy;
    M_All = sparse(NumNodes, NumFramesFinal);
    M_All(:,1) = M;
    F_All = sparse(NumNodes, NumFramesFinal);
    F_All(:,1) = F;
    Global_Area_All = sparse(NumNodes, NumFramesFinal);
    Global_Area_All(:,1) = Global_Area;
    Time_All = zeros(1, NumFramesFinal);
    Time_All(1,1) = SelectedTime;
    Area_All = zeros(1, NumFramesFinal);
    Area_All(1) = sum(Global_Area(InsideOldMask));
    Front_All = cell(1, NumFramesFinal);
    Front_All{1} = Front;
    FrontVelocity_All = cell(1, NumFramesFinal);
    FrontVelocity = zeros(size(Front.Polygon));
    for iVertex = 1:size(Front.Polygon,1)
        FrontVelocity(iVertex,:) = FT_velo_func(0, Front.Polygon(iVertex,:));
    end
    FrontVelocity_All{1} = FrontVelocity;
end

A_old = Global_Area(InsideOldMask);

fprintf('\n')
strlength = 0;
doVoronoi = false;
iNeigh = 1;

% segregated solver divergence detection
DivergenceCount = 0;

for iT = 1:NumTimeNodes
  
    if DivergenceCount > 10
      disp('Breaking the time iteration loop.');
      fprintf('Segregated solver diverged in the past %d iterations!',DivergenceCount);
       break;  
    end
    
    time = dt*iT;
    if ~mod(iT,round(NumTimeNodes/1000)) && ~config.DisplayFixedPointError
        strlength = DisplayProgress(ticID, iT, NumTimeNodes, strlength);
    end
    
    if ~config.FixedBoundary

        TotalArea = sum(A_old);
        if config.smoothing
            iNeigh = max(round(0.5*sqrt(size(Front.Polygon,1))),1);
        end

        Tangents = [Front.Polygon(end-iNeigh+1:end,:); Front.Polygon(1:end-iNeigh,:)] - [Front.Polygon(1+iNeigh:end,:); Front.Polygon(1:iNeigh,:)];
        Normals = [-Tangents(:,2), Tangents(:,1)]./repmat(sqrt(Tangents(:,1).^2+Tangents(:,2).^2),1,2);
        NormalX = scatteredInterpolant(Front.Polygon(:,1), Front.Polygon(:,2), Normals(:,1),'nearest');
        NormalY = scatteredInterpolant(Front.Polygon(:,1), Front.Polygon(:,2), Normals(:,2),'nearest');

        % different formulas for vp in ZS model
        if config.bilinearInterp
            Actin_Velocity = @(x,y) bilinear2D_vec(PosX, PosY, Ux, Uy, InteriorOldMask, x, y, MeshSize);
        else
            FX = scatteredInterpolant(PosX(InsideOldMask), PosY(InsideOldMask), Ux(InsideOldMask),'nearest');
            FY = scatteredInterpolant(PosX(InsideOldMask), PosY(InsideOldMask), Uy(InsideOldMask),'nearest');
        end        
        
    end
% Propogate front
    if ~config.FixedBoundary || iT==1
        [Front, FrontHitsOuterBox] = PropogateFront(f, time, Front, ExactPolygon, Domain, MeshSize);
        if FrontHitsOuterBox, break; end

    % Identify inside nodes    

        [InsideMask, InteriorMask, BoundaryMask] = IdentifyNodeTypes(PosX, PosY, Front, NumNodesX);

    % Update front movement metrics    
        [InteriorCandidatesMask, InteriorCandidatesIdx, BoundaryCandidatesIdx,new_nodes,gone_nodes] = ...
            FrontMovementUpdates(InsideMask, InteriorMask, BoundaryMask, InsideOldMask, InteriorOldMask, BoundaryOldMask);

    % Voronoi cell decomposition (Old nodes, new front)
    % Natural neighbors array: [i j], voronoi control volumes: A_i, edges: S_ij, inverse distances between nat.neighbors: D_Inv_ij 

        if doVoronoi || ~isempty(gone_nodes) || ~isempty(new_nodes)

            % update the discretization metrics that are dependent on voronoi
            % decomposition
            [NatNeighborBoundaryPairs, Vertices, VoronoiCells, D_Inv, Directions, MidPoints, NatNeighborPairs, Global_Area] = ...
                UpdateDiscretizationMetrics(BoundaryCandidatesIdx, PosX, PosY, InsideOldMask, NumNodesX, NumNodesY, InteriorCandidatesIdx, MeshSize, InteriorCandidatesMask);

        end
        % update the discretization metrics that are dependent on front
        % movement.
        [S, Global_Area, BoundaryIdxGlobal,BoundaryCurveLengths] = UpdateAreaEdgeLength(BoundaryCandidatesIdx, MeshSize, ...
            VoronoiCells, Vertices, NatNeighborBoundaryPairs, Directions, MidPoints, Front,true, true, InteriorCandidatesIdx, Global_Area);    


    % "Local" enumeration of the inside nodes 
        [I,J,NumInsideNodes] = LocalEnumeration(InsideOldMask,NumNodes,NatNeighborPairs);

    end
    
% Generate matrices for discretized system of equations       
% Update concentrations and velocities at the old nodes

    A = Global_Area(InsideOldMask);
    Ux_local = Ux(InsideOldMask);
    Uy_local = Uy(InsideOldMask);
    M_local = M(InsideOldMask);
    F_local = F(InsideOldMask);
    
    % --------------------------------------------------------------
    % preparing system matrices and RHS vectors for actin velocities
    % --------------------------------------------------------------
    % the viscous term should be assembled only once, if viscosity is
    % constant 
    if epsilon==0        % constant actin viscosity
        [Mat_Laplacian,~] = DiffusionTerm(eta, I, J, S, D_Inv, NumInsideNodes,true);
    end
   
    % the cell-substrate adhesion in actin velocity equations
    CellVariablekisi =  CreateCellVariable(ksi,NumInsideNodes);
    Mat_source = SourceTerm(A, CellVariablekisi); 
    
    % the active stress term in actin velocity equations
    [~,vector_x,~,vector_y] = GradientTerm(NatNeighborPairs, PosX, PosY, I, J, S, NumInsideNodes,sigma*M_local);
    
    % update right-hand-side vectors in actin velocity equations
    BUx = vector_x;
    BUy = vector_y;
    
    % the actin compressive stress (pressure) in actin velocity equations
    [~,vector_x,~,vector_y] = GradientTerm(NatNeighborPairs, PosX, PosY, I, J, S, NumInsideNodes,F_local);   
    % update right-hand-side vectors
    BUx = BUx - zeta*vector_x;
    BUy = BUy - zeta*vector_y;       
    
    % book-keeping for the old time step values
    M_local_PTS = M_local;
    F_local_PTS = F_local;
    
    % error printing initiation
    if config.DisplayFixedPointError
        if config.FixedPointStopCrite
            disp(['===================t=' num2str(SelectedTime+time) '==================='])
            fprintf('iter.\tUx\t\t\t     Uy\t\t\t    M\t\t\t    F\n');        
        end
    end

    
   % the segregated solver iteration loop
   for ss=1:config.FixedPointIterationsMax
       
      if config.FixedPointStopCrite
          % store the old values for error calculation
          M_local_old = M_local;
          F_local_old = F_local;
          Ux_local_old = Ux_local;
          Uy_local_old = Uy_local; 
      end
      
      % --------------------
      % solve for velocities
      % --------------------
      % if viscosity is non-constant it should be assembled at every
      % nonlinear solver iteration
      if epsilon ~=0                 % non-constant actin viscosity
        DiffCoeff = @(f) eta*(1 + epsilon*f/F0);
        [D_face] = DiffusionCoefficient(DiffCoeff,I,J,'HarmonicAverage',F_local);
        [Mat_Laplacian,~] = DiffusionTerm(D_face, I, J, S, D_Inv, NumInsideNodes,false);        
      end
      
      % the system matrix in actin velocity equations
      AU = Mat_Laplacian + Mat_source; 
      
      % solve the system
      Ux_local = AU\BUx;
      Uy_local = AU\BUy;
      
      
      % ------------------------------------------------
      % update advection velocities and solve for Myosin
      % ------------------------------------------------
      if config.AvoidSingbyAdvec
        detM_local = M_local<MaxMyosinAdvection;
        AdvectionField = {@(ux,uy,m) ux.*detM_local.*(1-m/MaxMyosinAdvection) @(ux,uy,m) uy.*detM_local.*(1-m/MaxMyosinAdvection)};
        EdgeVelocities = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference',Ux_local,Uy_local,M_local);
      else
        AdvectionField = {@(ux,uy) ux @(ux,uy) uy};
        EdgeVelocities = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference',Ux_local,Uy_local);
      end
      [Mat_advection,~] = AdvectionTerm(EdgeVelocities, I, J, S, NumInsideNodes); 
      
      
      % system matrix and right-hand-side vector of the myosin equation
      if config.AvoidSingbyDiffu   % non-linear diffusion coefficient
        detM_local = M_local<MaxMyosinDiffusion;
        DiffCoeff = @(m) D*detM_local.*(1-m/MaxMyosinDiffusion);
        [D_face] = DiffusionCoefficient(DiffCoeff,I,J,'HarmonicAverage',M_local);
        [Mat_diffusion,~] = DiffusionTerm(D_face, I, J, S, D_Inv, NumInsideNodes,false);

      else                         % constant diffusion coefficient
        
        [Mat_diffusion,~] = DiffusionTerm(D, I, J, S, D_Inv, NumInsideNodes,true);
          
      end
      
      if ss == 1
        [Mat_transient,BM] = TransientTerm(A, A_old, M_local_PTS, 'true', NumInsideNodes);
      end
      AM = Mat_transient + dt*(Mat_advection + Mat_diffusion);
      M_local = AM\BM;
      
      % ------------------------------------
      % update the Myosin contractile stress
      % ------------------------------------
      [~,vector_x,~,vector_y] = GradientTerm(NatNeighborPairs, PosX, PosY, I, J, S, NumInsideNodes,sigma*M_local);    
      
      % -------------------------------------------
      % update the right-hand-side velocity vectors
      % -------------------------------------------
      BUx = vector_x;
      BUy = vector_y;
      
      % ---------------
      % solve for actin
      % ---------------
      % advection term
      AdvectionField = {@(ux,uy) ux @(ux,uy) uy};
      EdgeVelocities = AdvectionVelocity(AdvectionField, NatNeighborPairs, PosX, PosY, I, J, NumInsideNodes, 'CentralDifference',Ux_local,Uy_local);
      [Mat_advection,~] = AdvectionTermUpwind(EdgeVelocities, I, J, S, NumInsideNodes); 
      
      % transient and the reaction terms    
      if strcmp(config.ActinBoundaryCond,'DirichletActin')
          if ss == 1
            [~,BF] = TransientTerm(A, A_old, F_local_PTS, 'false', NumInsideNodes);
            BoundaryIdxLocal = GetBoundaryIndexLocal(InsideOldMask,size(F_local),BoundaryIdxGlobal);
            RHS_F = ApplyDirichletBndCondVec(BF,F0,BoundaryIdxLocal);
          end
          % reaction matrix is the same as the transient term matrix, so no
          % need to build it again!  
          % system matrix
          AF = Mat_transient*(1+koff*dt) + dt*(Mat_advection);    
          [AF,~] = ApplyDirichletBndCondMat(AF,InsideOldMask,NumInsideNodes,size(F),BoundaryIdxGlobal);               
          
      else
          % correct for known flux from the membrane
          % the total flux is known: (F * u) * n | mem. = (F0 * u)*n
          if ss == 1
              [~,BF] = TransientTerm(A, A_old, F_local_PTS, 'false', NumInsideNodes);
              BoundaryIdxLocal = GetBoundaryIndexLocal(InsideOldMask,size(F_local),BoundaryIdxGlobal);
          end
          Flux = F0 *...
              (  NormalX(PosX(BoundaryIdxGlobal),PosY(BoundaryIdxGlobal)).*Ux_local(BoundaryIdxLocal) + ...
              NormalY(PosX(BoundaryIdxGlobal),PosY(BoundaryIdxGlobal)).*Uy_local(BoundaryIdxLocal)  ).* ...
              BoundaryCurveLengths;

          % add the flux to the RHS of advection equation
          RHS_F = BF;
          RHS_F(BoundaryIdxLocal) = RHS_F(BoundaryIdxLocal) - dt*Flux;
          % system matrix
          AF = Mat_transient*(1+koff*dt) + dt*(Mat_advection);
      end
      % solve the system
      F_local = AF\RHS_F;
      
      
      % ----------------------------------------------
      % update the actin compressive stress (pressure)
      % ----------------------------------------------
      [~,vector_x,~,vector_y] = GradientTerm(NatNeighborPairs, PosX, PosY, I, J, S, NumInsideNodes,F_local);   
      % update right-hand-side vectors
      BUx = BUx - zeta*vector_x;
      BUy = BUy - zeta*vector_y;      
      
      % Error control and stopping criteria
      if config.FixedPointStopCrite                % Maximum absolute error
            
          % calculate solution errors
          MaxErrorUx = full(max(abs(Ux_local-Ux_local_old)));
          MaxErrorUy = full(max(abs(Uy_local-Uy_local_old)));
          MaxErrorM = full(max(abs(M_local-M_local_old)));
          MaxErrorF = full(max(abs(F_local-F_local_old)));
          MaxErrorSol = max(max(max(MaxErrorUx,MaxErrorUy),MaxErrorM),MaxErrorF);
          if config.DisplayFixedPointError
              fprintf('%d\t%8.7e\t%8.7e\t%8.7e\t%8.7e\n',ss,MaxErrorUx,MaxErrorUy,MaxErrorM,MaxErrorF); 
          end          
          if MaxErrorSol<config.FixedPointIterationsTol
             break;
          end
          if ss==config.FixedPointIterationsMax && MaxErrorSol>=config.FixedPointIterationsTol
              fprintf('Warning: the segregated solver did not converge at t = %5.3f! Max. Sol. error is : %8.7e after %i iterations\n',SelectedTime+time,MaxErrorSol,config.FixedPointIterationsMax);
              DivergenceCount = DivergenceCount+1;
          end
          
      end

            
   end % end of segregated solver
    
     
    Ux(InsideOldMask) = Ux_local;
    Uy(InsideOldMask) = Uy_local;
    M(InsideOldMask)  = M_local;
    F(InsideOldMask)  = F_local;
   
% Adjust concentaration at the nodes neighboring those that disappear and new nodes
   if ~config.FixedBoundary
        [Ux, Uy, M, F, Global_Area, InsideOldMask, BoundaryOldMask, NatNeighborPairs] = GoneNodesAdjustmentActin(gone_nodes, Global_Area, Ux, Uy, M, F, PosX, PosY, InsideOldMask, BoundaryOldMask, NatNeighborPairs, NumNodesX, NumNodesY, MeshSize, Front);
        [Ux, Uy, M, F, Global_Area, InsideOldMask, NatNeighborPairs] = NewNodesAdjustmentActin(new_nodes, Global_Area, Ux, Uy, M, F, PosX, PosY, InsideOldMask, BoundaryOldMask, NatNeighborPairs, NumNodesX, NumNodesY, MeshSize, Front);
   end

% Collect the results    

    if config.returnAllData && ~mod(iT,max(round(NumTimeNodes/NumFrames),1))
        ii = iT/max(round(NumTimeNodes/NumFrames),1)+1;
        Time_All(ii) = SelectedTime+time;
        Area_All(ii) = TotalArea;
        Front_All{ii} = Front;
        FrontVelocity = zeros(size(Front.Polygon));
        for iVertex = 1:size(Front.Polygon,1)
            FrontVelocity(iVertex,:) = FT_velo_func(time, Front.Polygon(iVertex,:));
        end
        FrontVelocity_All{ii} = FrontVelocity;
        Ux_All(:,ii) = Ux;
        Uy_All(:,ii) = Uy;
        M_All(:,ii) = M;
        F_All(:,ii) = F;
        Global_Area_All(:,ii) = Global_Area;
        
        if config.SaveDuringSolve
          
          % Save into temporary data structures
          dataMask = Time_All>0;
          dataMask(1) = true;

          M_All_temp = M_All;
          F_All_temp = F_All;
          Ux_All_temp = Ux_All;
          Uy_All_temp = Uy_All;
          Time_All_temp = Time_All;
          Area_All_temp = Area_All;
          Front_All_temp = Front_All;
          FrontVelocity_All_temp = FrontVelocity_All;
          Global_Area_All_temp = Global_Area_All;

          M_All = M_All(:,dataMask);
          F_All = F_All(:,dataMask);
          Ux_All = Ux_All(:,dataMask);
          Uy_All = Uy_All(:,dataMask);
          Time_All = Time_All(dataMask);
          Area_All = Area_All(dataMask);
          Front_All = Front_All(dataMask);
          FrontVelocity_All = FrontVelocity_All(dataMask);
          Global_Area_All = Global_Area_All(:,dataMask);      

          save([CurrentDir delim_ch 'actin' delim_ch str '.mat'], 'config', 'M_tot', 'M_All', 'F_All', 'Ux_All', 'Uy_All', 'Time_All', 'Area_All', ...
          'PosX', 'PosY', 'Front_All', 'FrontVelocity_All',...
          'Global_Area_All', 'config','InitialMass','-v7.3', 'str');

          M_All = M_All_temp;
          F_All = F_All_temp;
          Ux_All = Ux_All_temp;
          Uy_All = Uy_All_temp;
          Time_All = Time_All_temp;
          Area_All = Area_All_temp;
          Front_All = Front_All_temp;
          FrontVelocity_All = FrontVelocity_All_temp;
          Global_Area_All = Global_Area_All_temp;

          clear M_All_temp F_All_temp Ux_All_temp Uy_All_temp Time_All_temp Area_All_temp Front_All_temp Global_Area_All_temp;
          
        end        
        
    end

    
    % modify for precise last time output    
    if config.returnAllData && (iT==iT_last_output) && det_last_output
        Time_All(ii+1) = SelectedTime+time;
        Area_All(ii+1) = TotalArea;
        Front_All{ii+1} = Front;
        FrontVelocity = zeros(size(Front.Polygon));
        for iVertex = 1:size(Front.Polygon,1)
            FrontVelocity(iVertex,:) = FT_velo_func(time, Front.Polygon(iVertex,:));
        end
        FrontVelocity_All{ii+1} = FrontVelocity;
        Ux_All(:,ii+1) = Ux;
        Uy_All(:,ii+1) = Uy;
        M_All(:,ii+1) = M;
        F_All(:,ii+1) = F;
        Global_Area_All(:,ii+1) = Global_Area;      
    end
    
% Visualize solution while running
    
    if config.Visualization && ~mod(iT,round(NumTimeNodes/100))
        if iT == round(NumTimeNodes/100)
            figure('Position', [100 400 1400 500]);
        end
        subplot(1,2,1)
        VisualizeSolution(F, PosX, PosY, MeshSize, SelectedTime+time, Front, 0, 1)
        axis equal
        xlim(Domain{1})
        ylim(Domain{2})
        colorbar
        subplot(1,2,2)
        VisualizeSolution(sqrt(Ux.^2+Uy.^2), PosX, PosY, MeshSize,SelectedTime+time, Front, 0, 0)
        colorbar
        hold on
        VisualizeSolution({Ux, Uy, 1}, PosX, PosY, MeshSize, SelectedTime+time, Front, 1,0,{'color','w','LineWidth',1})
        xlim(Domain{1})
        ylim(Domain{2})    
        hold off
        pause(0.1)
    end
    
    BoundaryOldMask = BoundaryMask;
    InteriorOldMask = InteriorMask;
    A_old = Global_Area(InsideMask);
    if isempty(gone_nodes) && isempty(new_nodes)
        doVoronoi = false;
    else
        doVoronoi = true;
    end
end;
fprintf('\n')

if ~config.returnAllData
    Ux_All = Ux;
    Uy_All = Uy;
    M_All = M;
    F_All = F;
    Time_All = SelectedTime+time;
    Area_All = TotalArea;
    Front_All{1} = Front;
    Global_Area_All = Global_Area;
    FrontVelocity_All = zeros(size(Front.Polygon));
    for iVertex = 1:size(Front.Polygon,1)
        FrontVelocity_All(iVertex,:) = FT_velo_func(time, Front.Polygon(iVertex,:));
    end
end

if isempty(ExactPolygon)
  % free internal front representation
  fcmd.command = 'free';
  fcmd.front = f;
  fronTierCmd(fcmd);
end

dataMask = Time_All>0;
dataMask(1) = true;
M_All = M_All(:,dataMask);
F_All = F_All(:,dataMask);
Ux_All = Ux_All(:,dataMask);
Uy_All = Uy_All(:,dataMask);
Time_All = Time_All(dataMask);
Area_All = Area_All(dataMask);
Front_All = Front_All(dataMask);
Global_Area_All = Global_Area_All(:,dataMask);

% save the output file
save([CurrentDir delim_ch 'actin' delim_ch str '.mat'], 'config', 'M_tot', 'M_All', 'F_All', 'Ux_All', 'Uy_All', 'Time_All', 'Area_All', 'PosX', 'PosY', 'Front_All', 'FrontVelocity_All', 'Global_Area_All','InitialMass', 'config','str', '-v7.3');
str = [CurrentDir delim_ch 'actin' delim_ch str ];
