Biharmonic equation with only zero normal derivative

Dear Freefem community,

I’m trying to manipulate the boundary conditions for the bilaplacian problem, like the code below.

load "Morley"
// Parameters
int nn = 10; real h = 0.01;
real f = 1;
// Mesh
mesh Th = square(nn, nn);
Th = adaptmesh(Th, h, IsMetric=1);
// Fespace
fespace Vh(Th, P2Morley); //The Morley finite element space 

Vh [u, ux, uy], [v, vx, vy];
// Macro
macro bilaplacien(u, v) (dxx(u)*dxx(v) + dyy(u)*dyy(v) + 2.*dxy(u)*dxy(v)) //
// Problem
solve bilap ([u, ux, uy], [v, vx, vy]) = int2d(Th)(
        bilaplacien(u, v)
    )
    - int2d(Th)(
        f*v
    )
    + on(1, 2, 3, 4, u=0, ux=0, uy=0)
    ;
// Plot
plot(u, cmm="u");

I’m wondering if is possible to remove the boundary condition u=0 and let only ux=uy=0. I tried on(1, 2, 3, 4, ux=0, uy=0) with no luck.

Any help is appreciated and thanks in advance!

An example is in the Freefem++ documentation, page 307 for Morley finite elements and you will see how the boundary conditions are taken into account https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/bilaplace-2d-PETSc.edp