// Implementation of the DDR method for the Kirchhoff--Love problem
// inspired by http://arxiv.org/abs/2507.17333

method ddr_kirchhoff_love_four_point_load {
  // Alternating four-point Gaussian loading (a smooth ``tablecloth'' load).
  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))
            )
    )

  // Homogeneous clamped boundary data.
  function zero_displacement(vector X) -> scalar = 0.0
  function zero_gradient(vector X) -> vector = vector(0.0, 0.0)

  // Discrete space
  space Uh {
    vertex Poly(0, scalar) called vertex_value
    vertex Poly(0, vector) called vertex_gradient
    edge Poly(k-1, scalar) called edge_value
    edge Poly(k, scalar) called edge_normal_derivative
    element Poly(k-1, scalar)
  }

  // Homogeneous clamped boundary conditions: u = 0 and grad(u) = 0.
  boundary conditions clamped_boundary_conditions on Uh {
    on vertex V:
      dof(V, vertex_value) = evaluate_at_vertex(zero_displacement)
    on vertex V:
      dof(V, vertex_gradient) = evaluate_at_vertex(zero_gradient)
    on edge E:
      dof(E, edge_value) = l2_project(zero_displacement, Poly(k-1, scalar))
    on edge E:
      dof(E, edge_normal_derivative) = l2_project(zero_displacement, Poly(k, scalar))
  }

  // Operators
  operator edge_hessian : Uh(u) -> Poly(k+1, vector) on edge E {
    forall w in Poly(k+1, vector):
      int(E) edge_hessian(u) dot w =
      int(E) dof(u, E, edge_value) * (tangential_derivative(tangential_derivative(w)) dot tangent(E))
      - int(E) dof(u, E, edge_normal_derivative) * (tangential_derivative(w) dot normal(E))
      + sum_vertices(orientation(V, E) * (
                                           dof(u, V, vertex_gradient) dot w(V)
                                           - dof(u, V, vertex_value) * (tangential_derivative(w) dot tangent(E))
                                         )
                    )
  }

  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, edge_value) * tangential_derivative(q)
        + sum_vertices(orientation(V, E) * dof(u, V, vertex_value) * q(V))
  }

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

  operator element_gradient: Uh(u) -> Poly(k, vector) on element T {
    forall w in Poly(k, vector):
      int(T) element_gradient(u) dot w =
      - int(T) dof(u, T) * div(w)
      + int(dT) edge_trace(u) * (w dot normal)
  }

  operator gradient_based_potential: Uh(u) -> Poly(k+1, scalar) on element T {
    forall w in CurlPolyComplement(k+2):
      int(T) gradient_based_potential(u) * div(w) =
      - int(T) element_gradient(u) dot w
      + int(dT) edge_trace(u) * (w dot normal)
  }

  operator edge_gradient : Uh(u) -> Poly(k, vector) on edge E {
    forall w in Poly(k, vector):
      int(E) edge_gradient(u) dot w =
      int(E) (
               edge_tangential_derivative(u) * (w dot tangent(E))
               + dof(u, E, edge_normal_derivative) * (w dot normal(E))
             )
  }

  operator edge_gradient_kp2_tangential_derivative: Uh(u) -> Poly(k+1, vector) on edge E {
    forall w in Poly(k+1, vector):
      int(E) edge_gradient_kp2_tangential_derivative(u) dot w =
      - int(E) edge_gradient(u) dot tangential_derivative(w)
      + sum_vertices(orientation(V, E) * dof(u, V, vertex_gradient) dot w(V))
  }

  operator edge_gradient_kp2: Uh(u) -> Poly(k+2, vector) on edge E {
    forall w in ZeroAveragePoly(k+3, vector):
      int(E) edge_gradient_kp2(u) dot tangential_derivative(w) =
      - int(E) edge_gradient_kp2_tangential_derivative(u) dot w
      + sum_vertices(orientation(V, E) * dof(u, V, vertex_gradient) dot w(V))
  }

  operator edge_trace_kp3 : Uh(u) -> Poly(k+3, scalar) on edge E {
    forall q in ZeroAveragePoly(k+4, scalar):
      int(E) edge_trace_kp3(u) * tangential_derivative(q) =
      - int(E) (edge_gradient_kp2(u) dot tangent(E)) * q
      + sum_vertices(orientation(V, E) * dof(u, V, vertex_value) * q(V))
  }

  operator element_hessian : Uh(u) -> Poly(k+1, matrix) on element T {
    forall zeta in Poly(k+1, matrix):
      int(T) element_hessian(u) dot zeta =
      int(T) dof(u, T) * div(div(zeta))
      - int(dT) edge_trace_kp3(u) * (div(zeta) dot normal)
      + int(dT) edge_gradient_kp2(u) dot (zeta dot normal)
  }

  operator hessian_based_potential: Uh(u) -> Poly(k+2, scalar) on element T {
    forall tau in CurlCurlPolyComplement(k+4):
      int(T) hessian_based_potential(u) * div(div(tau)) =
      int(T) element_hessian(u) dot tau
      - sum_edges(orientation(E, T) *
                     (
                       int(E) edge_gradient_kp2(u) dot (tau dot normal(E))
                       - int(E) edge_trace_kp3(u) * (div(tau) dot normal(E))
                     )
                 )
  }

  operator edge_difference : Uh(u) -> Poly(k+1, vector) on edge E of element T {
    forall v in Poly(k+1, vector):
      int(E) edge_difference(u) dot v =
      int(E) (element_hessian(u) dot tangent(E) - edge_hessian(u)) dot v
  }

  // Forms
  bilinear form kirchhoff_love : Uh(trial u) times Uh(test v) {
    sum_elements(
                  int(T) element_hessian(u) dot element_hessian(v)
                  + diameter(T) * int(dT) edge_difference(u) dot edge_difference(v)
                )
  }

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

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