using Gridap
model = DiscreteModelFromFile("../models/model.json")

writevtk(model,"model")

labels = get_face_labeling(model)

add_tag_from_tags!(labels,"diri0",["sides", "sides_c"])

add_tag_from_tags!(labels,"dirig",
  ["circle","circle_c", "triangle", "triangle_c", "square", "square_c"])

V0 = TestFESpace(
  reffe=:Lagrangian, order=1, valuetype=Float64,
  conformity=:H1, model=model, labels=labels,
  dirichlet_tags=["diri0", "dirig"])

g = 1
Ug = TrialFESpace(V0,[0,g])

using LinearAlgebra: norm
const p = 3
@law flux(∇u) = norm(∇u)^(p-2) * ∇u
f(x) = 1
res(u,v) = ∇(v)*flux(∇(u)) - v*f

@law dflux(∇du,∇u) =
  (p-2)*norm(∇u)^(p-4)*inner(∇u,∇du)*∇u + norm(∇u)^(p-2) * ∇du
jac(u,du,v) = ∇(v)*dflux(∇(du),∇(u))

trian = Triangulation(model)
degree=2
quad = CellQuadrature(trian,degree)
t_Ω = FETerm(res,jac,trian,quad)

op = FEOperator(Ug,V0,t_Ω)

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

import Random
Random.seed!(1234)
x = rand(Float64,num_free_dofs(Ug))
uh0 = FEFunction(Ug,x)
uh, = solve!(uh0,solver,op)

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