laplace-lagrange-mult.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

/*
   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);