using Gridap

u(x) = x[1] + x[2]
f(x) = 0

∇u(x) = VectorValue(1,1)

import Gridap: ∇
∇(::typeof(u)) = ∇u

∇(u) === ∇u

domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)

domain3d = (0,1,0,1,0,1)
partition3d = (4,4,4)
model3d = CartesianDiscreteModel(domain3d,partition3d)

writevtk(model,"model")

order = 1
V0 = TestFESpace(
  reffe=:Lagrangian, order=order, valuetype=Float64,
  conformity=:H1, model=model, dirichlet_tags="boundary")

U = TrialFESpace(V0,u)

trian = Triangulation(model)
degree = 2
quad = CellQuadrature(trian,degree)

a(u,v) = ∇(v)⊙∇(u)
b(v) = v*f

t_Ω = AffineFETerm(a,b,trian,quad)
op = AffineFEOperator(U,V0,t_Ω)

uh = solve(op)

e = u - uh

writevtk(trian,"error",cellfields=["e" => e])

l2(w) = w*w
h1(w) = a(w,w) + l2(w)

el2 = sqrt(sum( integrate(l2(e),trian,quad) ))

eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

tol = 1.e-10
@assert el2 < tol
@assert eh1 < tol

const k = 2*pi
u(x) = sin(k*x[1]) * x[2]
∇u(x) = VectorValue(k*cos(k*x[1])*x[2], sin(k*x[1]))
f(x) = (k^2)*sin(k*x[1])*x[2]

∇(::typeof(u)) = ∇u
b(v) = v*f

function run(n,order)

  domain = (0,1,0,1)
  partition = (n,n)
  model = CartesianDiscreteModel(domain,partition)

  V0 = TestFESpace(
    reffe=:Lagrangian, order=order, valuetype=Float64,
    conformity=:H1, model=model, dirichlet_tags="boundary")

  U = TrialFESpace(V0,u)

  trian = Triangulation(model)
  degree=2*order
  quad = CellQuadrature(trian,degree)

  t_Ω = AffineFETerm(a,b,trian,quad)
  op = AffineFEOperator(U,V0,t_Ω)

  uh = solve(op)

  e = u - uh

  el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
  eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

  (el2, eh1)

end

function conv_test(ns,order)

  el2s = Float64[]
  eh1s = Float64[]
  hs = Float64[]

  for n in ns

    el2, eh1 = run(n,order)
    h = 1.0/n

    push!(el2s,el2)
    push!(eh1s,eh1)
    push!(hs,h)

  end

  (el2s, eh1s, hs)

end

el2s, eh1s, hs = conv_test([8,16,32,64,128],2);

using Plots

plot(hs,[el2s eh1s],
    xaxis=:log, yaxis=:log,
    label=["L2" "H1"],
    shape=:auto,
    xlabel="h",ylabel="error norm")

function slope(hs,errors)
  x = log10.(hs)
  y = log10.(errors)
  linreg = hcat(fill!(similar(x), 1), x) \ y
  linreg[2]
end

slope(hs,el2s)

slope(hs,eh1s)