poly_rust

method hho_poisson {
  // Data corresponding to u(x,y) = sin(pi x) sin(pi y)
  function source(vector X) -> scalar =
    2.0 * pow(3.141592653589793, 2.0)
    * sin(3.141592653589793 * X[0]) * sin(3.141592653589793 * X[1])
  function boundary_data(vector X) -> scalar = 0.0

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

  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) dof(u, T) * div(grad(q))
      + int(dT) dof(u, E) * (grad(q) dot normal)
    constraint int(T) potential_reconstruction(u) = int(T) 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) - element_difference(u)) * q
  }

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

  bilinear form poisson : Uh(trial u) times Uh(test v) {
    sum_elements(
      int(T) grad(potential_reconstruction(u)) dot grad(potential_reconstruction(v))
      + pow(diameter(T), -1.0) * int(dT) edge_difference(u) * edge_difference(v)
    )
  }

  boundary conditions dirichlet_boundary_conditions on Uh {
    on edge E:
      dof(E) = l2_project(boundary_data, Poly(k, scalar))
  }

  linear problem hho_poisson_problem on Uh {
    lhs { poisson }
    rhs { load }
    boundary conditions dirichlet_boundary_conditions
    export { potential_reconstruction }
  }
}