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 45 46 47 48 49 50
|
load "Element_HCT"
load "qf11to25" // for tripleQF function
// Parameter
real f = 1;
// Mesh
mesh Th = square(40, 40); //mesh definition of Omega
// Fespaces
fespace Wh(Th, P2);
fespace Vh(Th, HCT); // HCT 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)) // end of macro
// Problem
// WARNING: the quadrature formula must be defined on 3 sub triangles
// the function tripleQF build this king of formula from classical quadrature
QF2 qfHCT5 = tripleQF(qf5pT);
solve bilap ([u, ux, uy], [v, vx, vy])
= int2d(Th, qft=qfHCT5)(bilaplacien(u, v))
- int2d(Th)(f*v)
+ on(1, 2, 3, 4, u=0, ux=0, uy=0)
;
// Plot
plot(u, cmm="u", wait=1, fill=1);
plot(ux, wait=1, cmm="u_x");
plot(uy, wait=1, cmm="u_y");
// Max & Error
Wh uu = u;
real umax = uu[].max;
int err = (abs(umax-0.0012782) > 1e-4);
cout << " uu max = " << umax << " ~ 0.0012782, err = " << err << endl;
// Plot
int n = 100, nn = n+10;
real[int] xx(nn), yy(nn);
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);
// End
assert(err == 0);
|