// a trick to build C1 finite element with P3 lagrangre finite element
// like HSIEH-CLOUCH-TOCHER finite element 
// the idee is insure the the jump (dn(u)) on all edges 
// by penalisatison ... 
// not to bad ... 
// F. Hecht juin 2014 ..

load "Element_P3" // for P3
load "splitmesh3" // to splite each triangle in 3 traingle

int n=100,nn=n+10;
real[int] xx(nn),yy(nn);

mesh Th=square(20,20);  // mesh definition of $\Omega$
//Th=splitmesh3(Th);
fespace Vh(Th,P3);      // finite element space

macro bilaplacien(u,v) ( dxx(u)*dxx(v)+dyy(u)*dyy(v)+2.*dxy(u)*dxy(v)) // fin macro 
real f=1;
Vh [u],[v];
real pena = 1e6;
macro dn(u) (dx(u)*N.x+dy(u)*N.y)//
solve bilap([u],[v]) =
    int2d(Th)(  bilaplacien(u,v) ) 
  + intalledges(Th) ( jump(dn(u))*jump(dn(v))*pena) 
   - int2d(Th)(f*v)
   
   + on(1,2,3,4,u=0) 
; 
   
plot(u,cmm="u", wait=1,fill=1);
real umax = u[].max; 
int err =  (abs(umax-0.0012782) > 1e-4); 
cout << " uu max " << umax << " ~  0.0012782  , err = " << err 
     << " " << abs(umax-0.0012782) <<endl; 


for (int i=0;i<=n;i++)
 {
   xx[i]=real(i)/n;
   yy[i]=u(0.5,real(i)/n); // value of uh at point (0.5, i/10.) 
 }
 plot([xx(0:n),yy(0:n)],wait=1);
 assert(err==0); 
