pkg load msh

fem_init_env ();

run lshape-domain.m;
mshd = fem_init_mesh (msho);

V = fem_fs_VelocityUpdate (mshd);
Q = fem_fs_PressureUpdate (mshd);

bcnoslip = fem_bc (V, @(x,y) [0;0], [3,4]);
bcout = fem_bc (Q, @(x,y) 0, 2);

t = 0;
dt = 0.01;
T = 3.;

k = fem_coeff ('k', @(x,y) dt);
f = fem_coeff ('f', @(x,y) [0; 0]);
u0 = fem_coeff ('u0', @(x,y) [0; 0]);
p0 = fem_coeff ('p1', @(x,y) 0);
A1 = fem_rhs_TentativeVelocity (V, k, bcnoslip);
A3 = fem_rhs_VelocityUpdate (V, bcnoslip);
# solve the problem for each time step
# We need to update only the lhs
while t < T
  t += dt;
  bcin = fem_bc (Q, @(x,y) sin(3.0*t), 1);

  b1 = fem_lhs_TentativeVelocity (V, k, u0, f, bcnoslip);

  utmp = A1 \ b1;
 
  u1 = fem_func ('u1', V, utmp);

  A2 = fem_rhs_PressureUpdate (Q, bcin, bcout);
  b2 = fem_lhs_PressureUpdate (Q, bcin, bcout, u1, k);

  ptmp = A2 \ b2;
  p1 = fem_func ('p1', Q, ptmp);

  b3 = fem_lhs_VelocityUpdate (V, k, u1, p1, bcnoslip);

  ut = A3 \ b3;
  u0 = fem_func ('u0', V, ut);

  #press Q to make the plot continue
  fem_plot (u0);
end