// Implementation of the HHO method for the Kirchhoff--Love problem
// from http://dx.doi.org/10.1051/m2an/2017065
// Visualization benchmark with an alternating four-point Gaussian load and homogeneous clamped boundary conditions.

method hho_kirchhoff_love_four_point_load {
  parameter pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286198

  // Alternating four-point Gaussian load (tablecloth test)
  parameter load_amplitude = 1000.0
  parameter load_width = 0.055
  parameter load_offset = 0.22

  function source(vector X) -> scalar = load_amplitude * (
      exp(-(
            pow(X[0] - (0.5 - load_offset), 2.0)
            + pow(X[1] - (0.5 - load_offset), 2.0)
          ) / (2.0 * pow(load_width, 2.0)))
    - exp(-(
            pow(X[0] - (0.5 + load_offset), 2.0)
            + pow(X[1] - (0.5 - load_offset), 2.0)
          ) / (2.0 * pow(load_width, 2.0)))
    - exp(-(
            pow(X[0] - (0.5 - load_offset), 2.0)
            + pow(X[1] - (0.5 + load_offset), 2.0)
          ) / (2.0 * pow(load_width, 2.0)))
    + exp(-(
            pow(X[0] - (0.5 + load_offset), 2.0)
            + pow(X[1] - (0.5 + load_offset), 2.0)
          ) / (2.0 * pow(load_width, 2.0)))
  )

  // Numerical parameters
  function stabilization_weight(T) -> scalar = n_element_edges(T) * pow(diameter(T), 2.0) * pow(element_area(T), -1.0)

  // Discrete space
  space Uh {
    edge Poly(k, scalar) called edge_value
    edge Poly(k, vector) called edge_gradient
    element Poly(k, scalar)
  }

  // Homogeneous clamped boundary conditions: u = 0 and grad(u) = 0
  function zero_scalar(vector X) -> scalar = 0.0
  function zero_vector(vector X) -> vector = vector(0.0, 0.0)

  boundary conditions clamped_boundary_conditions on Uh {
    on edge E:
      dof(E, edge_value) = l2_project(zero_scalar, Poly(k, scalar))
    on edge E:
      dof(E, edge_gradient) = l2_project(zero_vector, Poly(k, scalar))
  }

  // Operators
  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, edge_value) * (grad(q) dot normal)
    constraint int(T) potential_reconstruction(u) = int(T) dof(u, T)
  }

  operator deflection_reconstruction : Uh(u) -> Poly(k+2, scalar) on element T {
    forall w in orthogonal complement of Poly(1, scalar) relative to Poly(k+2, scalar):
      int(T) grad(grad(deflection_reconstruction(u))) dot grad(grad(w)) =
      int(T) dof(u, T) * div(div(grad(grad(w))))
      + int(dT) dof(u, E, edge_gradient) dot (grad(grad(w)) dot normal)
      - int(dT) dof(u, E, edge_value) * (div(grad(grad(w))) dot normal)

    forall w in Poly(1, scalar):
      int(T) deflection_reconstruction(u) * w =
      int(T) dof(u, T) * w
      // int(T) potential_reconstruction(u) * w
  }

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

  operator edge_gradient_difference : Uh(u) -> Poly(k, vector) on edge E of element T {
    forall tau in Poly(k, vector):
      int(E) edge_gradient_difference(u) dot tau =
      int(E) (grad(deflection_reconstruction(u)) - dof(u, E, edge_gradient)) dot tau
  }

  operator edge_value_difference : Uh(u) -> Poly(k, scalar) on edge E of element T {
    forall w in Poly(k, scalar):
      int(E) edge_value_difference(u) * w =
      int(E) (deflection_reconstruction(u) - dof(u, E, edge_value)) * w
  }


  // Forms
  bilinear form kirchhoff_love : Uh(trial u) times Uh(test v) {
    sum_elements(
                  int(T) grad(grad(deflection_reconstruction(u))) dot grad(grad(deflection_reconstruction(v)))
                  + stabilization_weight(T) * pow(diameter(T), -4.0) * int(T) element_difference(u) * element_difference(v)
                  + pow(diameter(T), -1.0) * int(dT) edge_gradient_difference(u) dot edge_gradient_difference(v)
                  + pow(diameter(T), -3.0) * int(dT) edge_value_difference(u) * edge_value_difference(v)
                )
  }

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

  // Linear problem
  linear problem hho_kirchhoff_love_four_point_load_problem on Uh {
    lhs {
      kirchhoff_love
    }
    rhs {
      load
    }
    boundary conditions clamped_boundary_conditions
    export { deflection_reconstruction }
  }
}