method hho_poisson {
  parameter pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198

  // Symbolic names for the integer boundary labels stored in the mesh
  boundary labels { left = 1, right = 2, bottom = 3, top = 4 }

  // Exact solution
  function exact_solution(vector X) -> scalar = sin(pi * X[0]) * sin(pi * X[1])

  function exact_solution_gradient(vector X) -> vector
  = vector(
            pi * cos(pi * X[0]) * sin(pi * X[1]),
            pi * sin(pi * X[0]) * cos(pi * X[1])
          )

  function exact_solution_flux(vector X) -> scalar on edge E =
    -exact_solution_gradient(X) dot normal(E)

  function source(vector X) -> scalar = 2.0 * pi * pi * exact_solution(X)

  // Functions for consistency tests
  function test_linear(vector X) -> scalar = X[0]- 4.0 * X[1] + 7.0
  function test_linear_gradient(vector X) -> vector = vector(1.0, -4.0)


  function test_quadratic(vector X) -> scalar = pow(X[0], 2.0) - 3.0 * X[0] * X[1] + pow(X[1], 2.0) - 5.0
  function test_quadratic_gradient(vector X) -> vector
  = vector(
            2.0 * X[0] - 3.0 * X[1],
            -3.0 * X[0] + 2.0 * X[1]
          )

  // Discrete space
  space Uh {
    element Poly(k, scalar)
    edge Poly(k, scalar)
  }

  // Interpolants
  interpolant interpolant_exact_solution on Uh {
    on element T:
      dof(T) = l2_project(exact_solution, Poly(k, scalar))
    on edge E:
      dof(E) = l2_project(exact_solution, Poly(k, scalar))
  }

  interpolant interpolant_test_linear on Uh {
    on element T:
      dof(T) = l2_project(test_linear, Poly(k, scalar))
    on edge E:
      dof(E) = l2_project(test_linear, Poly(k, scalar))
  }

  interpolant interpolant_test_quadratic on Uh {
    on element T:
      dof(T) = l2_project(test_quadratic, Poly(k, scalar))
    on edge E:
      dof(E) = l2_project(test_quadratic, Poly(k, scalar))
  }

  // Operators
  operator gradient_reconstruction : Uh(u) -> Poly(k, vector) on element T {
    forall tau in Poly(k, vector):
      int(T) gradient_reconstruction(u) dot tau =
      - int(T) dof(u, T) * div(tau)
      + sum_element_edges(int(E) dof(u, E) * (tau dot normal))

    test exactness for k = 0 against test_linear_gradient using interpolant_test_linear
    test exactness for k = 1 against test_quadratic_gradient using interpolant_test_quadratic
  }

  operator potential_reconstruction : Uh(u) -> Poly(k+1, scalar) on element T {
    forall q in Poly(k+1, scalar):
      int(T) grad(potential_reconstruction(u)) dot grad(q) =
      int(T) gradient_reconstruction(u) dot grad(q)
    constraint int(T) potential_reconstruction(u) = int(T) dof(u, T)

    test exactness for k = 0 against test_linear using interpolant_test_linear
    test exactness for k = 1 against test_quadratic using interpolant_test_quadratic
  }

  operator potential_field : Uh(u) -> Poly(k+1, scalar) on element T {
    potential_field(u) = dof(u, T)
  }

  operator element_difference : Uh(u) -> Poly(k, scalar) on element T {
    forall q in Poly(k, scalar):
      int(T) element_difference(u) * q =
      int(T) (potential_reconstruction(u) - dof(u, T)) * q
  }

  operator edge_difference : Uh(u) -> Poly(k, scalar) on edge E of element T {
    forall q in Poly(k, scalar):
      int(E) edge_difference(u) * q =
      int(E) (potential_reconstruction(u) - dof(u, E)) * q
  }

  // Forms
  bilinear form poisson : Uh(trial u) times Uh(test v) {
    sum_elements(
                  int(T) gradient_reconstruction(u) dot gradient_reconstruction(v)
                  + pow(diameter(T), -2.0) * int(T) element_difference(u) * element_difference(v)
                  + pow(diameter(T), -1.0)  * int(dT) edge_difference(u) * edge_difference(v)
                )
  }

  linear form load : Uh(test v) {
    sum_elements(int(T) source * dof(v, T))
  }

  linear form neumann_contribution : Uh(test v) {
    sum_boundary_edges(right)(int(E) exact_solution_flux * dof(v, E))
  }

  // Boundary conditions
  boundary conditions mixed_boundary_conditions on Uh {
    on edge E in left, bottom, top:
      dof(E) = l2_project(exact_solution, Poly(k, scalar))
  }

  // Norm functionals
  function h1_component_norm : Uh(v) -> scalar {
    sqrt(
          sum_elements(
                        int(T) squared_norm(grad(dof(v, T)))
                        + pow(diameter(T), -1.0) * int(dT) pow(dof(v, E) - dof(v, T), 2.0)
                      )
        )
  }

  function h1_reconstruction_norm : Uh(v) -> scalar {
    sqrt(poisson(v, v))
  }

  // Problem description
  linear problem hho_poisson_mixed_boundary_conditions_problem on Uh {
    lhs { poisson }
    rhs { load - neumann_contribution }

    boundary conditions mixed_boundary_conditions
    compute errors using interpolant_exact_solution {
      h1_component_norm,
      h1_reconstruction_norm
    }
    export { potential_reconstruction, potential_field }
  }
}