I decided to expand my earlier comment into an answer. I'd suggest using the Morley element that uses $P_2$ basis in each element and the degrees of freedom are
- values at vertices
- normal derivatives at edge midpoints
It is probably the easiest element for biharmonic problem implementation wise.
Here is an example using sp.fem library (version 8019aa7):
from spfem.mesh import MeshTri
from spfem.asm import AssemblerGlobal
from spfem.element import ElementGlobalMorley
from spfem.utils import direct,const_cell
import numpy as np
# define mesh
m=MeshTri()
m.refine(5)
# define element and assemble the stiffness matrix and load vector
e1=AbstractElementMorley()
a=AssemblerAbstract(m,e1)
def bilinear(ddu,ddv):
return ddu[0][0]*ddv[0][0]+\
ddu[0][1]*ddv[0][1]+\
ddu[1][0]*ddv[1][0]+\
ddu[1][1]*ddv[1][1]
def F(x):
return 50
A=a.iasm(bilinear)
f=a.iasm(lambda v,x: F(x)*v)
# find boundary and interior DOFs
bnd_nodes=m.boundary_nodes()
bnd_facets=m.boundary_facets()
bnd_facets_mid=0.5*(m.p[:,m.facets[0,bnd_facets]]+\
m.p[:,m.facets[1,bnd_facets]])
d1=a.dofnum_u.getdofs(N=bnd_nodes)
d2=a.dofnum_u.getdofs(F=bnd_facets)
I=np.setdiff1d(np.arange(a.dofnum_u.N),np.union1d(d1,d2))
# set boundary conditions
x=np.zeros(A.shape[0])
def G1(x,y):
return 0.05*np.sin(4*np.pi*y)
def G2(x,y):
return np.sin(4*np.pi*x)
# check orientation of boundary normals from sign of jacobian determinant
tris=np.ones(m.t.shape[1])
flip=np.nonzero((a.mapping.detDF(np.array([[0]]))<0).flatten())[0]
tris[flip]=-1.0
ori=tris[m.f2t[0,bnd_facets]]
x[d1]=G1(m.p[0,bnd_nodes],m.p[1,bnd_nodes])
x[d2]=ori*G2(bnd_facets_mid[0,:],bnd_facets_mid[1,:])
# solve the system
x=direct(A,f,I=I,x=x)
# draw the solution
xverts=x[a.dofnum_u.n_dof[0,:]]
m.plot3(xverts)
m.show()
The solution looks like this:

Note that both $g_1$ and $g_2$ are easily defined as Dirichlet-type boundary conditions.