% N+2 points along each axis
    N  = 10;
% elasticity
    r  = 1.0e-1;
% damping
    c  = 2.0e-1;
% downward force
    f0 = 1.0e-1;
% final time
    T  = 1000;
% number of time steps
    K  = 100000;
% time step
    ta = T/K;
% stopping threshold
    tol = 1.0e-4;

% repeatedly used values
    N1   = N+1;
    N2   = N+2;
    N1N1 = N1*N1;
    N1N2 = N1*N2;
    N2N2 = N2*N2;

% grid parameters
    h  = 1/N1;
    x  = linspace(0,1,N2)';
    y  = linspace(0,1,N2);

% gradient, directly from nodes to interfaces between them
    D1 = speye(N1);
    D2 = spdiags(kron([-1,1],ones(N2,1)),[0,1],N1,N2)/h;
    Gx = kron(D1,D2);
    D1 = speye(N2);
    D2 = spdiags(kron([-1,1],ones(N1,1)),[-1,0],N1,N1)/h;
    Gy = kron(D2,D1);
% gradient, including Dirichlet boundary
    D1 = speye(N2);
    D2 = spdiags(kron([-1,1],ones(N2,1)),[0,1],N1,N2)/h;
    Dx = kron(D1,D2);
    Dy = kron(D2,D1);

% gradient, averaging nodes to interfaces nearby, including Dirichlet boundary
    C1 = spdiags(kron([-1,1],ones(N2,1)),[-1,1],N2,N2)/(2*h);
    C1(1,2) = 1/h; C1(1,1) = -1/h; C1(N2,N2-1) = -1/h; C1(N2,N2  ) = 1/h;
    C2 = spdiags(kron([1,1],ones(N2,1)),[0,1],N1,N2)/2;
    Cx = kron(C2,C1);
    Cy = kron(C1,C2);

% downward force vector
    f1 = zeros(N2,N1);
    f2 = zeros(N2,N1);
    f3 = -f0*ones(N2,N1);
    f  = [f1(:);f2(:);f3(:)];

% initial state, u = (u0+u1,u2,u3), ut = (u1t,u2t,u3t)
    u0 = kron(x,ones(1,N2));
    u1 = zeros(N2,N2);       u1 = u1(:,2:N2); u1 = u1(:);
    u2 = kron(ones(N2,1),y); u2 = u2(:,2:N2); u2 = u2(:);
    u3 = zeros(N2,N2);       u3 = u3(:,2:N2); u3 = u3(:);

    u1t = zeros(N2,N2); u1t = u1t(:,2:N2); u1t = u1t(:);
    u2t = zeros(N2,N2); u2t = u2t(:,2:N2); u2t = u2t(:);
    u3t = zeros(N2,N2); u3t = u3t(:,2:N2); u3t = u3t(:);

    U   = [u1;u2;u3;u1t;u2t;u3t];

% start time stepping
    for k=1:K
        t = k*ta;

% impose fixed boundary
        u1 = U(1:N1N2);
            u1 = reshape(u1,N2,N1);
                u1 = [zeros(N2,1),u1] + u0;
                    u1 = u1(:);
        u2 = U((N1N2+1):(2*N1N2));
            u2 = reshape(u2,N2,N1);
                u2 = [zeros(N2,1),u2];
                    u2 = u2(:);
        u3 = U((2*N1N2+1):(3*N1N2));
            u3 = reshape(u3,N2,N1);
                u3 = [zeros(N2,1),u3];
                    u3 = u3(:);

% gradients on x interfaces
        u1x  = Dx*u1; u2x  = Dx*u2; u3x  = Dx*u3;
        u1y_ = Cy*u1; u2y_ = Cy*u2; u3y_ = Cy*u3;

        u1x = reshape(u1x,N1,N2); u1x = u1x(:,2:N2); u1x = u1x(:);
        u2x = reshape(u2x,N1,N2); u2x = u2x(:,2:N2); u2x = u2x(:);
        u3x = reshape(u3x,N1,N2); u3x = u3x(:,2:N2); u3x = u3x(:);

        u1y_ = reshape(u1y_,N1,N2); u1y_ = u1y_(:,2:N2); u1y_ = u1y_(:);
        u2y_ = reshape(u2y_,N1,N2); u2y_ = u2y_(:,2:N2); u2y_ = u2y_(:);
        u3y_ = reshape(u3y_,N1,N2); u3y_ = u3y_(:,2:N2); u3y_ = u3y_(:);

% gradients on y interfaces
        u1x_ = Cx*u1; u2x_ = Cx*u2; u3x_ = Cx*u3;
        u1y  = Dy*u1; u2y  = Dy*u2; u3y  = Dy*u3;

% area increments on x interfaces
        dEx = sqrt((u2x.*u3y_-u3x.*u2y_).^2 + ...
                   (u1x.*u3y_-u3x.*u1y_).^2 + ...
                   (u1x.*u2y_-u2x.*u1y_).^2);

% area increments on y interfaces
        dEy = sqrt((u2x_.*u3y-u3x_.*u2y).^2 + ...
                   (u1x_.*u3y-u3x_.*u1y).^2 + ...
                   (u1x_.*u2y-u2x_.*u1y).^2);

% coefficients on x and y interfaces
        cfx = r*(1-1./dEx);
        cfy = r*(1-1./dEy);

% coefficient matrix
%     utt + c ut = -Au + f
%     Au = r(1-1/|ux X uy|)[uy]'[uy]ux)_x
%        + r(1-1/|ux X uy|)[ux]'[ux]uy)_y
%     [u] = [  0,-u3, u2; u3,  0,-u1;-u2, u1,  0]
% [u]'[u] = [u2^2+u3^2,-u1u2,-u1u3;-u1u2,u1^2+u3^2,-u2u3;-u1u3,-u2u3,u1^2+u2^2]
        Uy =  [spdiags(u2y_.^2+u3y_.^2,0,N1N1,N1N1), ...
              -spdiags(u1y_.*u2y_,0,N1N1,N1N1), ...
              -spdiags(u1y_.*u3y_,0,N1N1,N1N1); ...
                -spdiags(u1y_.*u2y_,0,N1N1,N1N1), ...
                 spdiags(u1y_.^2+u3y_.^2,0,N1N1,N1N1), ...
                -spdiags(u2y_.*u3y_,0,N1N1,N1N1); ...
                  -spdiags(u1y_.*u3y_,0,N1N1,N1N1), ...
                  -spdiags(u2y_.*u3y_,0,N1N1,N1N1), ...
                   spdiags(u1y_.^2+u2y_.^2,0,N1N1,N1N1)];
        Ux =  [spdiags(u2x_.^2+u3x_.^2,0,N1N2,N1N2), ...
              -spdiags(u1x_.*u2x_,0,N1N2,N1N2), ...
              -spdiags(u1x_.*u3x_,0,N1N2,N1N2); ...
                -spdiags(u1x_.*u2x_,0,N1N2,N1N2), ...
                 spdiags(u1x_.^2+u3x_.^2,0,N1N2,N1N2), ...
                -spdiags(u2x_.*u3x_,0,N1N2,N1N2); ...
                  -spdiags(u1x_.*u3x_,0,N1N2,N1N2), ...
                  -spdiags(u2x_.*u3x_,0,N1N2,N1N2), ...
                   spdiags(u1x_.^2+u2x_.^2,0,N1N2,N1N2)];
        Ax  = kron(speye(3),Gx')* ...
              kron(speye(3),spdiags(cfx,0,N1N1,N1N1))*Uy ...
             *kron(speye(3),Gx);
        Ay  = kron(speye(3),Gy')* ...
              kron(speye(3),spdiags(cfy,0,N1N2,N1N2))*Ux ...
             *kron(speye(3),Gy);
        A   = Ax + Ay;

        B0  = [sparse(3*N1N2,3*N1N2),speye(3*N1N2);-A,-c*speye(3*N1N2)];

% update with Crank Nicholson
        B1  = speye(6*N1N2) - 0.5*ta*B0;
        B2  = speye(6*N1N2) + 0.5*ta*B0;

        U   = B1 \ (B2*U + ta*[zeros(3*N1N2,1);f]);
        U   = U.*(abs(U) > tol);

% graph the current state
        surf(reshape(u1,N2,N2), ...
             reshape(u2,N2,N2), ...
             reshape(u3,N2,N2));
        xlabel('u1'); ylabel('u2'); zlabel('u3');
%       view([135,30])
       view([30*t,30])
        e = (mean(dEx)+mean(dEy))/2;
        res = ta*norm(A*U(1:3*N1N2)-f)/norm(f);
        del = norm(U(3*N1N2+1:6*N1N2))/norm(U(1:3*N1N2));
        title(sprintf('time=%0.2f, area=%0.2f, res=%0.1e, del=%0.1e', ...
            t,e,res,del));
        grid on; axis([-0.6 1.6 -1.6 1.6 -3.1 0.1]); pbaspect([1 1 1]);
        drawnow;

        if ((res < tol) || (del < tol))
            break;
        end

    end
