1// a trick to build C1 finite element with P3 lagrangre finite element
2// like HSIEH-CLOUCH-TOCHER finite element
3// the idee is insure the the jump (dn(u)) on all edges
4// by penalisatison ...
5// not to bad ...
6// F. Hecht juin 2014 ..
7
8load "Element_P3" // for P3
9load "splitmesh3" // to splite each triangle in 3 traingle
10
11int n=100,nn=n+10;
12real[int] xx(nn),yy(nn);
13
14mesh Th=square(20,20);  // mesh definition of $\Omega$
15//Th=splitmesh3(Th);
16fespace Vh(Th,P3);      // finite element space
17
18macro bilaplacien(u,v) ( dxx(u)*dxx(v)+dyy(u)*dyy(v)+2.*dxy(u)*dxy(v)) // fin macro
19real f=1;
20Vh [u],[v];
21real pena = 1e6;
22macro dn(u) (dx(u)*N.x+dy(u)*N.y)//
23solve bilap([u],[v]) =
24    int2d(Th)(  bilaplacien(u,v) )
25  + intalledges(Th) ( jump(dn(u))*jump(dn(v))*pena)
26   - int2d(Th)(f*v)
27
28   + on(1,2,3,4,u=0)
29;
30
31plot(u,cmm="u", wait=1,fill=1);
32real umax = u[].max;
33int err =  (abs(umax-0.0012782) > 1e-4);
34cout << " uu max " << umax << " ~  0.0012782  , err = " << err
35     << " " << abs(umax-0.0012782) <<endl;
36
37
38for (int i=0;i<=n;i++)
39 {
40   xx[i]=real(i)/n;
41   yy[i]=u(0.5,real(i)/n); // value of uh at point (0.5, i/10.)
42 }
43 plot([xx(0:n),yy(0:n)],wait=1);
44 assert(err==0);