\begin{array}{rcl} -\Delta u & =& f \quad\mbox{on } \Omega,\\ \partial_n u&=& 0 \quad \mbox{on } \partial \Omega. \end{array}

\begin{array}{rcl} -\Delta u & =& f \quad\mbox{on } \Omega,\\ u&=& 0 \quad \mbox{on } \partial \Omega. \end{array}

verbosity=0;
//
// Mesh generation
//
int nsize=3;
mesh Th=square(nsize,nsize,[x*pi,y*pi]);
real hsize=pi/nsize;
fespace Vh(Th,P1);
//
// identity matrix of n=Vh.ndof size
//
int n = Vh.ndof;
real[int] one(n); one=1.;
matrix Id=one; //
set(Id,solver=sparsesolver);
macro grad(u) [dx(u),dy(u)]// 
//
// To generate finite element matrix with tgv at boundary dofs
//
varf a0(u,v) = int2d(Th)(grad(u)'*grad(v)) + on(1,2,3,4,u=0);
//
// To generate finite element matrix with no tgv at boundary dofs
// It amounts to a singular matrix because is the fem of
//
//       -Lap(u)=f con \partial_n u=0 on boundary
//
varf a(u,v) = int2d(Th)(grad(u)'*grad(v)) ;
//
// to generate a vector with 1 on boundary dofs
//
varf von(u,v)= on(1,2,3,4,u=1);
real[int] on1=von(0,Vh, tgv=1); // one on boundary dofs

int[int] I(n); // to get indices of boundary dofs

for (int k=0,i=0;i<n;++i)
 if( on1[i]) I[i]=-1; // set -1 at indices on boundary dofs
 else 
I[i]=k++; //numbers all remaining dofs starting at 0
//
// get matrix for problem with tgv at boundary dofs
//
matrix A0 = a0(Vh,Vh,solver=sparsesolver);
//
 // get matrix for singular problem (no tgv at boundary dofs)
 //
matrix A = a(Vh,Vh,solver=sparsesolver);
//
// Computational trick to get the cond(A) in L2 norm using the 
// eigenvalues
//
int nev=3; // the number of eigenvalues searched
real[int] evmin(nev),evmax(nev); 
//
// closest eig values to zero for A. Saved in evmin[0]
//
int k0=EigenValue(A,Id,sym=true,value=evmin);
//
// closest eig values to zero for A^-1. Saved in evmax[0]
//  As a consequence:
//
//                1/evmax[0] is the largest eig values for A
//
int k=EigenValue(Id,A,sym=true,value=evmax);
cout << " # of Converged eig values for A is = "<<k0 << " # of  Converged eig values for A^{-1} is =  " << k << endl;
cout << 1/evmin[0] << " " << evmax[0] << endl;
cout << "  Pure Neumann problem (singular) " <<endl;
//
// L^2 condition number=  (largest eig value)/(min eig value)=(1/evmax[0]) / evmin[0]
//
cout << "  cond number of A = " << 1/(evmax[0] * evmin[0]) << endl;
//
//  Remove from A the rows and columns for boundary dofs
//  i.e., with tgv=1 on diagonal entry
//
cout<<"---------------------------------------------------------------------- "<<endl;
cout<<"Vector I sets -1 at boundary dofs and numbers the rest starting from 0 "<<endl;
cout<<"---------------------------------------------------------------------- "<<endl;
cout<<I<<endl;
matrix A00 = A(I^-1,I^-1);
//
//  Remove from A0 the rows and columns for boundary dofs, 
//  i.e., with tgv=1e30 on diagonal entry 
//  It serves to check result
//
matrix A01 = A0(I^-1,I^-1);
matrix Id0 = Id(I^-1,I^-1); 
//
// set solvers
//
set(A00,solver=sparsesolver);
set(A01,solver=sparsesolver);
set(Id0,solver=sparsesolver);

//same trick as before for A00

k0=EigenValue(A00,Id0,sym=true,value=evmin);
k =EigenValue(Id0,A00,sym=true,value=evmax);
cout << " # of Converged eig values for A00 is = "<<k0 << " # of Converged eig values for A00^{-1} is =  " << k << endl;
cout << k0 << " " << k << endl;
cout << 1/evmin[0] << " " << evmax[0] << endl;
cout << "  Homogeneous Dirichlet problem (no singular) " <<endl;
cout << "  Matrix deduced from that of pure Neumann problem " <<endl;
real trucond= 1/(evmax[0] * evmin[0]);
cout << " True cond number  = " << trucond << endl;
cout << " h  = " <<hsize << endl;
cout<<" cond= O(h^-2) ==> cond*h^2 approx constant "<<endl;
cout << " True cond number*h^2  = " << trucond*hsize^2 << endl;

//same trick as before for A01

k0=EigenValue(A01,Id0,sym=true,value=evmin);
k =EigenValue(Id0,A01,sym=true,value=evmax);
cout << "# of Converged eig values for A01 is = "<<k0 << "# of Converged eig values for A01^{-1} is =  " << k << endl;
cout << k0 << " " << k << endl;
cout << 1/evmin[0] << " " << evmax[0] << endl;
cout << "  Homogeneous Dirichlet problem (no singular) " <<endl;
cout << "  Matrix generated from that of Dirichlet problem " <<endl;
cout << " True cond number  = " << 1/(evmax[0] * evmin[0]) << endl;