function coupling
% COUPLING OF AN ELASTIC ELLIPSE WITH A FLUID. 
% IMPLICIT ON DIFFUSION. CONSTANT TIME STEP.
%
%    Uses level-set formulation of immersed boundary method from
%    Cottet & Maitre, M3AS, 2006.
%
% 03/2008 by Emmanuel Maitre
% Feel free to modify for teaching and learning.
%
%    The WENO scheme for advection equation on level set is taken from 
%    Baris Sumengen level-set package
%    http://vision.ece.ucsb.edu/~sumengen/
%    Two functions from this package have been inserted at the end
%    of this file, to provide a self contained code.
%    
%    The fluid solver from B. Seibold solves the incompressible 
%    Navier-Stokes equations in a rectangular domain with 
%    prescribed velocities along the boundary.
%    The solution method is finite differencing on
%    a staggered grid with implicit diffusion and a Chorin
%    projection method for the pressure.
%    Visualization is done by a colormap-isoline plot for
%    pressure and normalized quiver and streamline plot for
%    the velocity field.
%
% 07/2007 by Benjamin Seibold
% Feel free to modify for teaching and learning.
%-----------------------------------------------------------------------
Re = 1e2;     % Reynolds number
dt = 1e-2;    % time step WARNING THIS MUST BE CONSTANT
              % TO PERFORM CHOLESKI FACTORISATION OUT 
              % OF TIME LOOP.
tf = 2.2e0;    % final time
lx = 2;       % width of box
ly = 2;       % height of box
nx = 64;      % number of x-gridpoints
ny = 64;      % number of y-gridpoints
nsteps = 5;  % number of steps with graphic output
A = 0.75;     % horizontal radius of ellipse (fraction of lx)
B = 0.5;      % vertical radius of ellipse (fraction of ly)
lambda = 1;   % stiffness of elastic ellipse.
stretch = 1.26253110; % Initial stretching of elastic ellipse.
% The elastic response is encoded in functions Ep and its derivative Epp 
% Default is linear response to strain (see end of file).

%-----------------------------------------------------------------------
nt = ceil(tf/dt); dt = tf/nt;
x = linspace(0,lx,nx+1); hx = lx/nx;
y = linspace(0,ly,ny+1); hy = ly/ny;
eps = 1.5*hx;
%-----------------------------------------------------------------------
% initial conditions
% Velocity
U = zeros(nx-1,ny); V = zeros(nx,ny-1); 

fprintf('initialization')
% Initial curve in level-set
phi = zeros(nx,ny);
for i=1:nx
  for j = 1:ny
      phi(i,j)=d((x(i)+x(i+1))/2,(y(j)+y(j+1))/2,A,B);
  end
end

% Initialisation at initial stretch. Stretching is recorded in |grad phi|
% See Cottet & Maitre (M3AS, 2006)
phi = stretch * phi; 

% boundary conditions
uN = x*0;    vN = avg(x)*0;
uS = x*0;      vS = avg(x)*0;
uW = avg(y)*0; vW = y*0;
uE = avg(y)*0; vE = y*0;
%-----------------------------------------------------------------------
Ubc = dt/Re*([2*uS(2:end-1)' zeros(nx-1,ny-2) 2*uN(2:end-1)']/hx^2+...
      [uW;zeros(nx-3,ny);uE]/hy^2);
Vbc = dt/Re*([vS' zeros(nx,ny-3) vN']/hx^2+...
      [2*vW(2:end-1);zeros(nx-2,ny-1);2*vE(2:end-1)]/hy^2);

Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye(nx));
Lp(1,1) = 3/2*Lp(1,1);
perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp';
Lu = speye((nx-1)*ny)+dt/Re*(kron(speye(ny),K1(nx-1,hx,2))+...
     kron(K1(ny,hy,3),speye(nx-1)));
peru = symamd(Lu); Ru = chol(Lu(peru,peru)); Rut = Ru';
Lv = speye(nx*(ny-1))+dt/Re*(kron(speye(ny-1),K1(nx,hx,3))+...
     kron(K1(ny-1,hy,2),speye(nx)));
perv = symamd(Lv); Rv = chol(Lv(perv,perv)); Rvt = Rv';
Lq = kron(speye(ny-1),K1(nx-1,hx,2))+kron(K1(ny-1,hy,2),speye(nx-1));

perq = symamd(Lq); Rq = chol(Lq(perq,perq)); Rqt = Rq';

fprintf(', time loop\n--20%%--40%%--60%%--80%%-100%%\n')
for k = 1:nt
    % treat nonlinear terms
   gamma = min(1.2*dt*max(max(max(abs(U)))/hx,max(max(abs(V)))/hy),1);
   Ue = [uW;U;uE]; Ue = [2*uS'-Ue(:,1) Ue 2*uN'-Ue(:,end)];
   Ve = [vS' V vN']; Ve = [2*vW-Ve(1,:);Ve;2*vE-Ve(end,:)];
   Ua = avg(Ue')'; Ud = diff(Ue')'/2;
   Va = avg(Ve);   Vd = diff(Ve)/2;
   UVx = diff(Ua.*Va-gamma*abs(Ua).*Vd)/hx;
   UVy = diff((Ua.*Va-gamma*Ud.*abs(Va))')'/hy;
   Ua = avg(Ue(:,2:end-1));   Ud = diff(Ue(:,2:end-1))/2;
   Va = avg(Ve(2:end-1,:)')'; Vd = diff(Ve(2:end-1,:)')'/2;
   U2x = diff(Ua.^2-gamma*abs(Ua).*Ud)/hx;
   V2y = diff((Va.^2-gamma*abs(Va).*Vd)')'/hy;
   U = U-dt*(UVy(2:end-1,:)+U2x);
   V = V-dt*(UVx(:,2:end-1)+V2y);
   
    % compute curvature
    phixx = diff(phi,2)/hx^2; phixx = phixx([1 1:end end],:);
    phiyy = diff(phi',2)'/hy^2; phiyy = phiyy(:,[1 1:end end]);
    phix = (phi(3:end,:)-phi(1:end-2,:))/(2*hx); phix = phix([1 1:end end],:);
    phiy = (phi(:,3:end)-phi(:,1:end-2))/(2*hy); phiy = phiy(:,[1 1:end end]);
    phixy = (phix(:,3:end)-phix(:,1:end-2))/(2*hy); phixy = phixy(:,[1 1:end end]);
    grad = phix.^2+phiy.^2;
    courb = (phixx.*phiy.^2-2*phix.*phiy.*phixy+phiyy.*phix.^2)./grad.^1.5;
    prot = (-(phixx.*phix+phixy.*phiy).*phiy+(phixy.*phix+phiyy.*phiy).*phix)./grad; % VALID ONLY IF STRETCH = ngp    


   % compute elastic force
    dxphi = diff(phi)/hx; dyphi = diff(phi')'/hy;
    dyphi = dyphi(1:end-1,:)+dyphi(2:end,:);
    dyphi = dyphi(:,1:end-1)+dyphi(:,2:end);
    dyphi = dyphi/4;
   
    ngp = sqrt(dxphi(:,2:end-1).^2 + dyphi.^2); ngp=ngp(:,[1 1:end end]);
   
    F = (prot(1:end-1,:)+prot(2:end,:))/2 .* Epp(lambda,ngp) .* (-dyphi(:,[1 1:end end]))./ngp; 
    F = F - Ep(lambda,ngp).*(courb(1:end-1,:)+courb(2:end,:))/2 .* dxphi;     
    
    U = U + dt * F./(eps*ngp) .* zeta((phi(1:end-1,:)+phi(2:end,:))/2/eps./ngp);
   
    dyphi = diff(phi')'/hy; dxphi = diff(phi)/hx;
    dxphi = dxphi(:,1:end-1)+dxphi(:,2:end);
    dxphi = dxphi(1:end-1,:)+dxphi(2:end,:);
    dxphi = dxphi/4;
   
    ngp = sqrt(dyphi(2:end-1,:).^2 + dxphi.^2); ngp=ngp([1 1:end end],:);
    
    F = (prot(:,1:end-1)+prot(:,2:end))/2 .* Epp(lambda,ngp) .* dxphi([1 1:end end],:)./ngp; 
    F = F - Ep(lambda,ngp).*(courb(:,1:end-1)+courb(:,2:end))/2 .* dyphi;  
   
    V = V + dt * F./(eps*ngp) .* zeta((phi(:,1:end-1)+phi(:,2:end))/2/eps./ngp);
  
   % implicit viscosity
   rhs = reshape(U+Ubc,[],1);
   u(peru) = Ru\(Rut\rhs(peru));
   U = reshape(u,nx-1,ny);
   rhs = reshape(V+Vbc,[],1);
   v(perv) = Rv\(Rvt\rhs(perv));
   V = reshape(v,nx,ny-1);
   
   % pressure correction
   rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1);
   p(perp) = -Rp\(Rpt\rhs(perp))/dt;
   P = reshape(p,nx,ny);
   U = U-dt*diff(P)/hx;
   V = V-dt*diff(P')'/hy;
   
   % move curve (transport equation)
    phi = phi - dt * evolve_vector_WENO(phi,hy,hx,[vS' (V(:,1:end-1)+V(:,2:end))/2 vN'],[uE; (U(1:end-1,:)+U(2:end,:))/2; uW]);   
   
   % visualization
   if floor(25*k/nt)>floor(25*(k-1)/nt), fprintf('.'), end
   if k==1|floor(nsteps*k/nt)>floor(nsteps*(k-1)/nt)
      figure(1)
      % stream function
      rhs = reshape(diff(U')'/hy-diff(V)/hx,[],1);
      q(perq) = Rq\(Rqt\rhs(perq));
      Q = zeros(nx+1,ny+1);
      Q(2:end-1,2:end-1) = reshape(q,nx-1,ny-1);
      % pressure
      clf, contourf(avg(x),avg(y),P',50,'LineColor', 'none'), hold on
      contour(x,y,Q',20,'w-'),
      colormap jet
      % zero level set
      contour(avg(x),avg(y),phi',[0 0],'k-');
      % velocity field
      %Ue = [uS' avg([uW;U;uE]')' uN'];
      %Ve = [vW;avg([vS' V vN']);vE];
      %Len = sqrt(Ue.^2+Ve.^2+eps);
      %quiver(x,y,(Ue./Len)',(Ve./Len)',.4,'k-')
      hold off, axis equal, axis([0 lx 0 ly])
      p = sort(p); caxis(p([8 end-7])); colorbar;
      title(sprintf('Pressure and stream function at time t = %0.2g',k*dt))
      figure(2)
      plot(avg(x),(P(:,ny/2)+P(:,ny/2+1))/2-(P(1,ny/2)+P(1,ny/2+1))/2,'-o');
      title(sprintf('Pressure profile at time t = %0.2g',k*dt))
      drawnow
   end
end

fprintf('\n')

%====================== FUNCTIONS FOR NS SOLVER ===========================

function B = avg(A,k)
if nargin<2, k = 1; end
if size(A,1)==1, A = A'; end
if k<2, B = (A(2:end,:)+A(1:end-1,:))/2; else B = avg(A,k-1); end
if size(A,2)==1, B = B'; end

function A = K1(n,h,a11)
% a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
A = spdiags([-1 a11 0;ones(n-2,1)*[-1 2 -1];0 a11 -1],-1:1,n,n)'/h^2;

%====================== FUNCTIONS FOR LEVEL SET ===========================

function z = zeta(r) % cut-off function
   z = (r>-1).*(r<1).*(0.5*(1+cos(pi*r)));

     
% Compute signed distance to initial curve (only for initialization)
function res = d(xi,yj,A,B)
   N = 2400; ds = 2*pi/N; s = (0:N)*ds;   
   [dmin,is] = min((xs(s,A)-xi).^2+(ys(s,B)-yj).^2);
   tx = (xs(s(is+1),A)-xs(s(is-1),A))/(2*ds);
   ty = (ys(s(is+1),B)-ys(s(is-1),B))/(2*ds);
   if ((xi-xs(s(is),A))*ty - (yj-ys(s(is),B))*tx>0) 
       res = sqrt(dmin); else res = -sqrt(dmin);
   end
   
% The following test case is from Lee & Leveque IIM method (2002)
function res = xs(s,A)
   res = 1+A*cos(s);
   
function res = ys(s,B)
   res = 1+B*sin(s);
        
        
function res = Ep(lmbd, str) % Elastic response law
   % Interface fluide
   %res = lmbd*ones(size(str));
   % Interface elsatique 
   res = lmbd * (str-1);
function res = Epp(lmbd, str) % Derivative of Elastic response law wrt its second argument
   % Interface fluide
   %res = 0;
   % Interface elsatique 
   res = lmbd;

   
% ========================== WENO ADVECTION ===============================

function [delta] = evolve_vector_WENO(phi, dx, dy, u_ext, v_ext)
%
% Finds the amount of evolution under a vector field
% based force and using 5th order accurate WENO scheme
%
% Author: Baris Sumengen  sumengen@ece.ucsb.edu
% http://vision.ece.ucsb.edu/~sumengen/
%

delta = zeros(size(phi)+6);
data_ext = zeros(size(phi)+6);
data_ext(4:end-3,4:end-3) = phi;
vit = zeros(size(phi)+6);
vit(4:end-3,4:end-3) = u_ext;
% first scan the rows
for i=1:size(phi,1)
    delta(i+3,:) = delta(i+3,:) + upwind_WENO(data_ext(i+3,:), vit(i+3,:), dx);
end
% then scan the columns
vit(4:end-3,4:end-3) = v_ext;
for j=1:size(phi,2)
    delta(:,j+3) = delta(:,j+3) + upwind_WENO(data_ext(:,j+3), vit(:,j+3), dy);
end
delta = delta(4:end-3,4:end-3);

function [data_x] = upwind_WENO(data, F, dx)
%
% Calculates the upwind derivative using
% fifth order accurate WENO scheme
% takes 1-D data
% data: input data
% F: force field
% dx: grid resolution
% Note: before entering this function, data and F need to be 
% extended by 3 at the beginning and end (values don't matter)
%
% Author: Baris Sumengen  sumengen@ece.ucsb.edu
% http://vision.ece.ucsb.edu/~sumengen/
%


if length(F) == 1
    F = F*ones(size(data));
end

if length(data) ~= length(F)
    error('force and data lengths does not match');
end

data_x = zeros(size(data));

% extrapolate the beginning and end points of data
data(3) = 2*data(4)-data(5);
data(2) = 2*data(3)-data(4);
data(1) = 2*data(2)-data(3);
data(end-2) = 2*data(end-3)-data(end-4);
data(end-1) = 2*data(end-2)-data(end-3);
data(end) = 2*data(end-1)-data(end-2);

D1 = (data(2:end)-data(1:end-1))/dx;

for i=1:(length(data)-6)
    if F(i+3) > 0  % use D-
        k = i-1;
        v1 = D1(k+1);
        v2 = D1(k+2);
        v3 = D1(k+3);
        v4 = D1(k+4);
        v5 = D1(k+5);
    elseif F(i+3) < 0  % use D+
        k = i;
        v1 = D1(k+5);
        v2 = D1(k+4);
        v3 = D1(k+3);
        v4 = D1(k+2);
        v5 = D1(k+1);
    else
        continue;
    end

    data_x_1 = v1/3 - 7*v2/6 + 11*v3/6;
    data_x_2 = -v2/6 + 5*v3/6 + v4/3;
    data_x_3 = v3/3 + 5*v4/6 - v5/6;

    epsilon = 1e-6*max([v1 v2 v3 v4 v5].^2) ...
        + 1e-99;

    S1 = (13/12)*(v1-2*v2+v3).^2 ...
        + 0.25*(v1-4*v2+3*v3).^2;
    S2 = (13/12)*(v2-2*v3+v4).^2 ...
        + 0.25*(v2-v4).^2;
    S3 = (13/12)*(v3-2*v4+v5).^2 ...
        + 0.25*(3*v3-4*v4+v5).^2;

    alpha1 = 0.1/((S1+epsilon).^2);
    alpha2 = 0.6/((S2+epsilon).^2);
    alpha3 = 0.3/((S3+epsilon).^2);
    a_total = alpha1+alpha2+alpha3;
    data_x(i+3) = (alpha1/a_total)*data_x_1 + (alpha2/a_total)*data_x_2 ...
        + (alpha3/a_total)*data_x_3;

    data_x(i+3) = F(i+3)*data_x(i+3);
end


