load "Element_P3"
load "MUMPS_seq"

//defaulttoMUMPSseq();
verbosity=4;
 mesh Th=square(10,10);
 fespace Vh(Th,P3);     // P1 FE space
 Vh  uh,vh;              // unkown and test function.
 func f=1;                 //  right hand side function
 func g=0;                 //  boundary condition function
 real[int] rinfo(40);
 int[int] info(40);
 varf laplace(uh,vh) =                    //  Definion of  the problem
    int2d(Th,qft=qf7pT)( dx(uh)*dx(vh) + dy(uh)*dy(vh) - 100*uh*vh   ) //  bilinear form
  + int2d(Th,qft=qf7pT)( f*vh )                          //  linear form
  + on(1,2,3,4,uh=g) ;                      //  boundary condition form
matrix A=laplace(Vh,Vh,solver=sparsesolverSym);
set(A,solver=sparsesolver,rinfo=rinfo,info=info);
real[int] b=laplace(0,Vh);
cout << " nnz = "<< A.nnz << endl;  
//verbosity=1;
  uh[]=A^-1*b; // solve the problem plot(uh); // to see the result
  cout << "uh[] " << uh[].linfty << endl;
  plot(uh,ps="o/LaplaceP3.eps",value=true);
  if(verbosity>1)
 for(int i=0; i<40;++i)
   {
    cout << " solve: info " << i+1 << " = " << info[i] << " , / r " << rinfo[i] <<endl;
   }
