pkg load msh
pkg load fem-fenics

fem_init_env ();

fem_create_all ("TentativeVelocity");
fem_create_all ("VelocityUpdate");
fem_create_all ("PressureUpdate");

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

V = VelocityUpdate_FunctionSpace (mshd);
Q = PressureUpdate_FunctionSpace (mshd);

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

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

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

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

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

  A2 = PressureUpdate_BilinearForm (Q, bcin, bcout);
  b2 = PressureUpdate_LinearForm (Q, bcin, bcout, u1, k);

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

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

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

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