% Solution of the nonlinear problem
%
% \nabla \cdot (Vp F)+s0*F=0, 								(1)
% \nabla \cdot Vf = 0,        								(2)
% \ksi0 F (Vf - Vp) + \nabla p = 0, 						(3)
% \nabla(Vf*a - D\nabla a)-k*(b-a)-s0*F = 0,                (4)
% \nabla(Vf*a - D\nabla b)+k*(b-a) = 0,                     (5)
%
%
% where \nabla={d/dx, d/dy}
% Vp={-x/(1+exp(3*y)), -1/(1+exp(-3*y))} = {Vpx, Vpy}.
% Vf={Vfx, Vfy} is 2D vector fluid velocity field,
% F is filament(F)-actin density,
% a is the G-actin-profilin density, and
% b is the G-actin-thymosin density,
% p is fluid pressure.
%
%
% The boundary conditions are
%  Variable       LS          RS          FE            RE
%     F           F0          F0         f(x)        no-flux
%     Vf                 NO FLUX     \nabla \cdot Vf = 0
%     p                       FREE 
%    fa        -F*(Vp)                -F*(Vp)            0
%     fb                  NO FLUX     = 0
%
%   Here (fa) = Vf*a - D\nabla a, (fb) = Vf*b - D\nabla b,
%   F0=0.5, f(x)=1-x*x/8, k=5, D =1.
%
%
% A straightforward solution is not possible due to bad coditioned
% matrix, so we need to use several tricks to reach the solution.
%
% First, the F-actin problem is solved separately
% in the structure aF. This solution is then used as input data
% for the water velocity problem. To do this we try advanced approach
% using the function femval_stat.m 

% The problem to be solved is
% \nabla \cdot Vf = 0,        							(2)
% F*(Vf - Vp) + 1/ksi0 \nabla p = 0, 					(3)
%
% and it is presented in the coefficient form.
%
% 
% The aF mode is written in the coefficient form with flux term, i.e., 
% \nabla \cdot (-alpha FF) + a FF = f, where 
% a=s0, f=-eps*(s0+Vpxx+Vpyy), alpha={-Vpx -Vpy}.
% The Dirichlet BC read:
% LS and RS - FF=F0-eps;
% FE        - FF=f(x)-eps.
% No-flux BC on RE is (nx*FF_x+ny*FF_y) = 0,
% which is coded using
% q=-(nx*Vpx+ny*Vpy), and g=(nx*FF_x+ny*FF_y).

% The no-flux BC should be specified in the form
% nx*u+ny*v = 0 only for velocity components {u,v}. This implies
% g={u*dnx+v*dny, u*dnx+v*dny, 0} at all boundaries.

% Constants    s0 = 1, D=1, \ksi0 = 0.3.
%
% Solution of the hydrodynamics problem is performed in the 
% another auxiliary structure aFw, which uses the F-actin density
% data obtained at the previous step.

% When we are done with preparatory work we may start the solution
% of the reaction-diffusion-convection problem for G-actin specified
% in the structure fem. We again try to access data from aF and aFw using the
% functions 'femval_stat.m' for the quantities defined on the subdomains, and
% 'femval_statBnd.m' which works for the boundary-defined variables.
%
% To visualize the field Vp use the following commands  
% [xx,yy]=meshgrid(-2:0.2:2,-1:0.1:1);
% vvpx=-xx./(1+exp(3*yy)); vvpy=-1./(1+exp(-3*yy));
% This one shows the arrows for the vector field.
% quiver(xx,yy,vvpx,vvpy);
% hold on;
%
%    

clear fem aF aFw
% Spatial variables
fem.sdim={'x' 'y'};
% fem.shape=2;
% Constants and domain data
th0=pi/8;
a=2.6;
b=0.5;
hmin=0.06;  % size of the triangle
arscale = 1.3; % arrow scale for the water velocity plot
xscale = 1.0;
% domain geometry
clear s c p

%%
%% Domain definition
%%
 % upper ellipse
 ellipse1=ellip2(0,0,a,a,0);
 % lower ellipse
 ellipse2=ellip2(0,-0.0*b,a,b,0);
 rec=rect2(-a,a,-a,0);
 ro=ellipse1-rec-ellipse2;
 
 % ellipse2=ellip2(0,-2.0*b+1.6*b,0.55*a,0.50*b,0); %deep pit
 % ellarge=ellip2(0,-2.0*b+1.2*b,a,0.3*b,0);
 % ro=ellipse1-ellipse2-ellarge;
 
 % This removes all internal boundaries
 r=geomdel(ro);
 % r=ellip2(0,0,1.5,0.5,0);
% Define objects
objs={r};
names={'r'};
s.objs=objs;
s.name=names;
% 
% Curve objects
%
objs={};
names={};
c.objs=objs;
c.name=names;

objs={};
names={};
p.objs=objs;
p.name=names;

drawstruct=struct('s',s,'c',c,'p',p);
fem.draw=drawstruct;
fem.geom=geomcsg(fem);
% This enables to see the numbering of the boundary 
% segments of the domain.
%%%%%
geomplot(fem, 'edgelabels', 'on', 'edgearrow', 'on');axis equal;
pause;
%%%%%
% It appears that the following indices are given:
%
%
% RE = 1 & 3
% FE = 2 & 4
%
% Specify the F-actin mode 
aF.sdim={'x' 'y'};
aF.dim={'F'};
aF.shape=2;
% Boundary conditions are Dirichlet on all segments
% There are three types of boundary segments:
% RS and LS are of the type 1 (Dirichlet nonzero constant),
% FE is of the type 3 (Dirichlet nonzero variable) 
%
% aF.bnd.ind = [1 1 1 1 1 1];
aF.bnd.ind = [1 2 1 2];
% Constants
aF.const={'V0x',0.0,'V0y',-0.90, ...
        'Vpmin',0.1,'Vpmax',1.0,'Fmax',0.8,'Fmin',0.20, ...
        'aa',a,'bb',b,'tt',th0};

% The coefficient form
aF.form='coefficient';
aF.equ.a='s0';
aF.equ.f='0';
aF.equ.al={{{'-Vpx' '-Vpy'}}};
% aF.bnd.h=1;
% aF.bnd.r='prof';
aF.bnd.h={1 0};
aF.bnd.r={'prof' 0};
aF.bnd.q={'0' '-(dnx*Vpx+dny*Vpy)'};
aF.bnd.g={'0' '(dnx*Fx+dny*Fy)'};

aF.geom=fem.geom;
% Make initial mesh
aF.mesh=meshinit(aF,'hmax',hmin);
aF.expr = {'s0' 2.8 'F0' 1 ...
        'Vpx' 'V0x-x*vel0(x,y,aa,bb,Vpmax,Vpmin,tt)/sqrt(x*x+y*y)' ...
        'Vpy' 'V0y-y*vel0(x,y,aa,bb,Vpmax,Vpmin,tt)/sqrt(x*x+y*y)' ...
        'prof' 'profile(x,y,aa,bb,Fmax,Fmin)'};

aF.equ.init={{'F0'}};
% aF.equ.da = {{1}};
aF.xmesh = meshextend(aF); 
% Solve the problem for F-actin.

aF.sol = femlin(aF);

% Plot F-actin density.
postplot(aF,'tridata','prof','tribar','on',...
'arrowdata',{'Vpx' 'Vpy'},'arrowcolor','black',...
'Title','F-actin velocity');
axis equal;
pause;
% postplot(aF,'tridata','Vpy','tribar','on','Title','y vel');
%axis equal;
%pause;
postplot(aF,'tridata','F','tribar','on',...
    'Title','F-actin density');
daspect([1 xscale 1]); pbaspect([xscale 1 1]);
% axis equal;
% axes('PlotBoxAspectRatio',[1 2 1]);
pause;

% Define aFw mode (F-actin + water)
aFw.sdim={'x' 'y'};
aFw.dim = {'u','v','p'};
aFw.expr = {'ksi0' 0.3 'eps' 0.000000};
% aFw.const={'eps',0.0};
aFw.geom=fem.geom;
% Make initial mesh
aFw.mesh=meshinit(aFw,'hmax',hmin);

aFw.bnd.ind = [1 1 1 1 1 1];
aFw.shape = {shlag(1,'p') shlag(2,'u')  shlag(2,'v')};

aFw.form = 'coefficient';
aFw.equ.a = {{'femval_stat(x,y,''F'',''aF'')';...
            'femval_stat(x,y,''F'',''aF'')'; 0}};  % F*Vf
aFw.equ.be = {{ ...
{0 0} {0 0} {'1/ksi0/femval_stat(x,y,''F*F'',''aF'')' 0}; ...
{0 0} {0 0} {0 '1/ksi0/femval_stat(x,y,''F*F'',''aF'')'}; ...
{0 0} {0 0} {0 0} ...
}};   %   1/ksi0 \nabla p  = \beta \cdot \nabla u
aFw.equ.al = {{ ...
 {0 0} {0 0} {0 0}; ...
 {0 0} {0 0} {0 0}; ...
 {-1 0} {0 -1} {0 0} ...
}};    % \nabla \cdot Vf  = -\nabla \cdot \alpha u
aFw.equ.f={{...
             'femval_stat(x,y,''F*Vpx'',''aF'')' ...
             'femval_stat(x,y,''F*Vpy'',''aF'')' ...
             '0' ...
         }};

% No-flux BC
 
aFw.bnd.q = {{...
 'dnx' 'dny' '0'; ...
 'dnx' 'dny' '0'; ...
 '0' '0' '0' ...
}};
%  aFw.bnd.g = {{'0' ... 'u*dnx+v*dny' ...
%                '0' ... 'u*dnx+v*dny' ...
%                '0'}};

% Zero velocity BC
% aFw.bnd.r={{0 0 0}};
% aFw.bnd.h={{1 1 1}};


% For the general form
% aFw.form='general';
% aFw.equ.ga={{ ...
%            {0 0} ... 
%            {0 0} ...
%            {'u' 'v'} ...
%        }};
% aFw.equ.f={{ ...
%            '-(u-Vpx)*F-px/ksi0' ...
%            '-(v-Vpy)*F-py/ksi0' ...
%            '0' ...
%        }};
% aFw.bnd.g={ ...
%            {'u*dnx+v*dny' 'u*dnx+v*dny' '0'} ... % LS and RS
%            {'u*dnx+v*dny' 'u*dnx+v*dny' '0'} ... % FE
%            {'u*dnx+v*dny' 'u*dnx+v*dny' '0'} ... % RE
%        };


% aFw=femdiff(aFw);   % For general form only
aFw.xmesh = meshextend(aFw); 
% Solve the equations and plot the result.
aFw.sol = femlin(aFw);

% Plot the pressure and show the water velocity field on the same 
% picture.
postplot(aFw,'tridata','p','tribar','on',...
   'Title','pressure + water velocity',...
   'arrowdata',{'u' 'v'},'arrowcolor','black','arrowscale',arscale); 
%'arrowstyle','normalized');
daspect([1 xscale 1]); pbaspect([xscale 1 1]);
% axis equal;
pause;
% postplot(aFw,'tridata','ux+vy','tribar','on',...
%   'Title','velocity divergence');
% axis equal;
% pause;
% Plot separately the pressure gradient field.
% postarrow(fem,{'px' 'py'},'Title','Pressure gradient','Arrowcolor','r');
% or pressure gradient with the water velocity fields
% postplot(fem,'arrowdata',{'px' 'py'},...
%   'Title','pressure gradient + water velocity',...
%   'arrowdata',{'u' 'v'},'arrowcolor','black');
% axis equal;
%pause;
% Checking the accuracy of the solution - plot the 
% scalar norm of the LHS of eq. (3).
% postplot(aFw,'tridata',...
%'sqrt((-(u-Vpx)*femval_stat(x,y,''F'',''aF'')-px/ksi0)^2+(-(v-Vpy)*femval_stat(x,y,''F'',''aF'')-py/ksi0)^2)',...
%'tribar','on',...
%'Title','Error',...
%'Tridlim',[0 0.1]);
% axis equal;
% pause;

% Define the fem mode for solution of the G-actin problem
fem.sdim={'x' 'y'};
% Constants and domain data
fem.expr = {'D' 1 'k' 5 'A0' 0.3 'B0' 0.1};
fem.dim = {'A','B'};
fem.bnd.ind = [1 1 1 1 1 1];
fem.shape = 2;
% The coefficient form
fem.form='coefficient';

%
% The trick is not to use the original form of the equation with
% the source term s0*F. Instead of it we replace this term by the
% divergence of some flux term, namely s0*F=-\nabla \cdot (Vp*F).
%
% This trick leads to resolving of the positive source term problem
% which otherwise would lead to the solution growing in time.




% fem.equ.f = {{'femval_stat(x,y,''F*s0'',''aF'')' '0'}};  % {s0*F, 0} in RHS
fem.equ.f = {{'0' '0'}};  % {s0*F, 0} in RHS
fem.equ.c = {{'D' '0'; '0' 'D'}}; % -D \nabla^2 {a ,b} -> 
%  \nabla \cdot (-c \nabla U)
fem.equ.al = {{ ...
{'femval_stat(x,y,''-u'',''aFw'')' 'femval_stat(x,y,''-v'',''aFw'')'} {0 0}; ...  \% \nabla \cdot (Vf*{a, b})
{0 0} {'femval_stat(x,y,''-u'',''aFw'')' 'femval_stat(x,y,''-v'',''aFw'')'}...    % \nabla \cdot (-alpha u)
}};
fem.equ.ga = {{ ...
{'femval_stat(x,y,''Vpx*F'',''aF'')' 'femval_stat(x,y,''Vpy*F'',''aF'')'}; ...
{'0' '0'} ...
}};
fem.equ.a={{'k' '-k'; '-k' 'k'}};  % {k(a-b), -k(a-b)} -> a U
%
% BC is \n \cdot (c \nabla U + alpha U - gamma) = g
%
% fem.bnd.q={...
%    '0'; ... 
%    '0'; ... 
%    '0' ... 
% };

% fem.bnd.g={ ... 
%       {'0' '0'} ... % {'femval_statBnd(dom,s,''F*(Vpx*dnx+Vpy*dny)'',''aF'')' '0'} ...  % LS and RS
%       {'0' '0'} ... % {'femval_statBnd(dom,s,''F*(Vpx*dnx+Vpy*dny)'',''aF'')' '0'} ...  % FE        
%       {'0' '0'} ... % RE %  
%  };

fem.equ.init={{'A0' 'B0'}};
% fem.equ.da = 0.1;

fem.mesh=meshinit(fem,'hmax',hmin);
%fem=femdiff(fem);   % For general form only
fem.xmesh = meshextend(fem); 
% Solve the equations and plot the result.



fem.sol = femlin(fem);
% fem.sol=femtime(fem,'tlist',linspace(0,2.0,12),'report','on');
postplot(fem,'tridata','A','tribar','on',...
   'Title','G-actin-profilin'); 
daspect([1 xscale 1]); pbaspect([xscale 1 1]);
% axis equal;
pause;
postplot(fem,'tridata','B','tribar','on',...
   'Title','G-actin-thymosin'); 
daspect([1 xscale 1]); pbaspect([2 xscale 1]);
% axis equal;
pause;
% postplot(fem,'tridata',... 
%    'femval_stat(x,y,''ux+vy'',''aFw'')*A +(femval_stat(x,y,''u'',''aFw'')*Ax+femval_stat(x,y,''v'',''aFw'')*Ay)-D*(Axx+Ayy)-femval_stat(x,y,''s0*F'',''aF'')-k*(A-B)',...
%    'tribar','on',...
%   'Title','Error G-profilin'); axis equal;
% pause;
% postplot(fem,'tridata',... 
%    'femval_stat(x,y,''ux+vy'',''aFw'')*B +(femval_stat(x,y,''u'',''aFw'')*Bx+femval_stat(x,y,''v'',''aFw'')*By)-D*(Bxx+Byy)+k*(A-B)',...
%    'tribar','on',...
%   'Title','Error G-thymosin'); axis equal;
% pause;
postplot(fem,...
    'tridata',... 
    'A+B+femval_stat(x,y,''F'',''aF'')',...
    'tribar','on',...
   'Title','Total actin'); 
daspect([1 xscale 1]); pbaspect([xscale 1 1]);
% axis equal;
pause;
% plot of the profilin density on the front edge
postcrossplot(fem, 1, [2 5], 'lindata','A','cont','on','lincolor','r');
hold on;
% plot of the thymosin density on the front edge
postcrossplot(fem, 1, [2 5],'lindata','B','cont','on','lincolor','b',...
   'Title','G-actin');
hold off;