1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
|
// 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);
|
