code: def local_project(v, V, u=None):
      
      dv = TrialFunction(V)
      v_ = TestFunction(V)
      a_proj = inner(dv, v_)*dx
      b_proj = inner(v, v_)*dx
      solver = LocalSolver(a_proj, b_proj)
      solver.factorize()
      if u is None:
          u = Function(V)
          solver.solve_local_rhs(u)
          return u
      else:
          solver.solve_local_rhs(u)
          return
docs: Element-wise projection using LocalSolver