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
|
/*
solving laplace operator with neumam boundary condition
with 1D lagrange multiplier
The variationnal form is
find (u,l) such that
$\forall (v,m) a(u,v) + b(u,m) + b(v,l) = L(v) $
where $b(u,m) = int u*m dx$
*/
mesh Th=square(10,10);
fespace Vh(Th,P1); // P1 FE space
int n = Vh.ndof;
int n1 = n+1;
Vh uh,vh; // unkown and test function.
func f=1+x-y; // right hand side function
varf va(uh,vh) = // definion of the problem
int2d(Th)( dx(uh)*dx(vh) + dy(uh)*dy(vh) ) // bilinear form
;
varf vL(uh,vh)= int2d(Th)( f*vh ) ;
varf vb(uh,vh)= int2d(Th)(1.*vh);
matrix A=va(Vh,Vh);
real[int] b(n);
b = vL(0,Vh);
real[int] B = vb(0,Vh);
// the block matrix
matrix AA = [ [ A , B ] ,
[ B', 0 ] ] ;
real[int] bb(n+1),x(n+1),b1(1),l(1);
b1=0;
// build the block rhs
bb = [ b, b1];
set(AA,solver=UMFPACK);
x = AA^-1*bb; // solve the linear systeme
[uh[],l] = x; // set the value
cout << " l = " << l(0) << " , b(u,1) =" << B'*uh[] << endl;
plot(uh,wait=1);
|
