// Implementation of a variation of the DDR method for the Poisson problem
// from http://arxiv.org/abs/2605.23405
// The variation lies in the stabilization bilinear form

method ddr_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)

  // Discrete space
  space Uh {
    vertex Poly(0, scalar)
    edge Poly(k-1, scalar)
    element Poly(k-1, scalar)
  }

  // Interpolant
  interpolant interpolant_exact_solution on Uh {
    on vertex V:
      dof(V) = evaluate_at_vertex(exact_solution)
    on edge E:
      dof(E) = l2_project(exact_solution, Poly(k-1, scalar))
    on element T:
      dof(T) = l2_project(exact_solution, Poly(k-1, scalar))
  }

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

  // Operators
  operator edge_tangential_derivative : Uh(u) -> Poly(k, scalar) on edge E {
    forall q in Poly(k, scalar):
      int(E) edge_tangential_derivative(u) * q =
        - int(E) dof(u, E) * tangential_derivative(q)
        + sum_vertices(orientation(V, E) * dof(u, V) * q(V))
  }

  operator edge_potential_reconstruction : Uh(u) -> Poly(k+1, scalar) on edge E {
    forall q in ZeroAveragePoly(k+2, scalar):
      int(E) edge_potential_reconstruction(u) * tangential_derivative(q) =
        - int(E) edge_tangential_derivative(u) * q
        + sum_vertices(orientation(V, E) * dof(u, V) * q(V))
  }

  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)
        + int(dT) edge_potential_reconstruction(u) * (tau dot normal)
  }

  operator potential_reconstruction : Uh(u) -> Poly(k+1, scalar) on element T {
    forall tau in CurlPolyComplement(k+2):
      int(T) potential_reconstruction(u) * div(tau) =
        - int(T) gradient_reconstruction(u) dot tau
        + int(dT) edge_potential_reconstruction(u) * (tau dot normal)
  }

  operator edge_difference : Uh(u) -> Poly(k+1, scalar) on edge E of element T {
    forall q in Poly(k+1, scalar):
      int(E) edge_difference(u) * q =
        int(E) (potential_reconstruction(u) - edge_potential_reconstruction(u)) * 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), -1) * sum_element_edges(int(E) edge_difference(u) * edge_difference(v))
    )
  }

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

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

  // Error norms functionals
  function h1_potential_norm : Uh(v) -> scalar {
    sqrt(poisson(v, v))
  }

  function l2_potential_norm : Uh(v) -> scalar {
    sqrt(
          sum_elements(
                        int(T) pow(potential_reconstruction(v), 2.0)
                        + diameter(T) * int(dT) pow(edge_difference(v), 2.0)
                      )
        )
  }

  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, T) - dof(v, E), 2.0)
                        + sum_vertices(pow(dof(v, T) - dof(v, V), 2.0))
                      )
        )
  }

  function l2_component_norm : Uh(v) -> scalar {
    sqrt(
          sum_elements(
                        int(T) pow(dof(v, T), 2.0)
                        + diameter(T) * int(dT) pow(dof(v, E), 2.0)
                        + pow(diameter(T), 2.0) * sum_vertices(pow(dof(v, V), 2.0))
                      )
        )
  }

  // Problem description
  linear problem ddr_poisson_problem on Uh {
    lhs {
      poisson
    }
    rhs {
      load - neumann_contribution
    }
    boundary conditions mixed_boundary_conditions
    export { potential_reconstruction }
    compute errors using interpolant_exact_solution {
      h1_potential_norm, l2_potential_norm, h1_component_norm, l2_component_norm
    }
  }
}