File: bilapHCT.edp

package info (click to toggle)
freefem%2B%2B 3.61.1%2Bdfsg1-4
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 17,108 kB
  • sloc: cpp: 141,214; ansic: 28,664; sh: 4,925; makefile: 3,142; fortran: 1,171; perl: 844; awk: 290; php: 199; pascal: 41; f90: 32
file content (50 lines) | stat: -rw-r--r-- 1,188 bytes parent folder | download | duplicates (3)
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);