using Gridap
n = 100
domain = (0,1,0,1)
partition = (n,n)
model = CartesianDiscreteModel(domain,partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"diri1",[6,])
add_tag_from_tags!(labels,"diri0",[1,2,3,4,5,7,8])

D = 2
order = 2
V = TestFESpace(
  reffe=:Lagrangian, conformity=:H1, valuetype=VectorValue{D,Float64},
  model=model, labels=labels, order=order, dirichlet_tags=["diri0","diri1"])

Q = TestFESpace(
  reffe=:PLagrangian, conformity=:L2, valuetype=Float64,
  model=model, order=order-1, constraint=:zeromean)

uD0 = VectorValue(0,0)
uD1 = VectorValue(1,0)
U = TrialFESpace(V,[uD0,uD1])
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V, Q])
X = MultiFieldFESpace([U, P])

const Re = 10.0
@law conv(u,∇u) = Re*(∇u')⋅u
@law dconv(du,∇du,u,∇u) = conv(u,∇du)+conv(du,∇u)

function a(x,y)
  u, p = x
  v, q = y
  ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u)
end

c(u,v) = v⊙conv(u,∇(u))
dc(u,du,v) = v⊙dconv(du,∇(du),u,∇(u))

function res(x,y)
  u, p = x
  v, q = y
  a(x,y) + c(u,v)
end

function jac(x,dx,y)
  u, p = x
  v, q = y
  du, dp = dx
  a(dx,y)+ dc(u,du,v)
end

trian = Triangulation(model)
degree = (order-1)*2
quad = CellQuadrature(trian,degree)
t_Ω = FETerm(res,jac,trian,quad)
op = FEOperator(X,Y,t_Ω)

using LineSearches: BackTracking
nls = NLSolver(
  show_trace=true, method=:newton, linesearch=BackTracking())
solver = FESolver(nls)

uh, ph = solve(solver,op)

writevtk(trian,"ins-results",cellfields=["uh"=>uh,"ph"=>ph])