// different ways to solve the same problem in the same way border Gamma(t = 0,2*pi) { x = cos(t); y = sin(t); } mesh Th = buildmesh(Gamma(300)); fespace Vh(Th,P1); int n = Vh.ndof; real[int] phiv(n),rhsv(n); Vh phi,w, f = 1; // Standard Laplace, with CG solver solve Laplace(phi,w,solver = CG) = int2d(Th)(dx(phi)*dx(w)+dy(phi)*dy(w)) -int2d(Th)(f*w)+on(Gamma,phi = 1); cout << phi[].l2 << endl; // set the problem and later solve problem Laplace1(phi,w,solver = CG) = int2d(Th)(dx(phi)*dx(w)+dy(phi)*dy(w)) -int2d(Th)(f*w)+on(Gamma,phi = 1); Laplace1; cout << phi[].l2 << endl; // now I try to solve explicitely the linear system // bilinear form varf vLaplace(phi,w) = int2d(Th)(dx(phi)*dx(w)+dy(phi)*dy(w)) +on(Gamma,phi = 1); // compute once and for all the associated matrix matrix A = vLaplace(Vh,Vh,solver=CG); // specify the solver // the right hand side: it is a bilinear form with a dummy argument varf rhs(unused,w) = int2d(Th)(f*w)+on(Gamma,unused = 1); rhsv = rhs(0,Vh); // it is possbile, for small numbers, to display the matrix and the rhs // cout << A << endl; // cout << rhsv << endl; phiv = A^(-1)*rhsv; // no inversion of A cout << phiv.l2 << endl; // now LinearCG // now the func that applies the matrix to a vector func real[int] Laplace2phi(real[int] & phi) { real[int] r = A*phi; return r; } // diagonal preconditioner real[int] prec(n); // extract the diagonal of A prec = A.diag; // the func that applies the preconditioner func real[int] P(real[int] & phi) { phi ./= prec; return phi; } // initial vector for LinearCG (it has to satisfy bc) // trick to force bc varf boundary(u,v) = on(1,u=1); real[int] onboundary = boundary(0,Vh); phiv = 0; // enforce bc phiv = onboundary ? 1. : phiv; // try to comment this LinearCG(Laplace2phi,phiv,rhsv,nbiter=1000,precon=P); cout << phiv.l2 << endl; // now I try without the costruction of the matrix. Of course, it is slower func real[int] newLaplace2phi(real[int] & phi) { Vh philoc; philoc[] = phi; varf newvLaplace(unused,w) = int2d(Th)(dx(philoc)*dx(w)+dy(philoc)*dy(w)) +on(Gamma,unused = 1); // you can use or not bc // if you omit bc, then you NEED the preconditioner real[int] r = newvLaplace(0,Vh); return r; } phiv = 0; // enforce bc phiv = onboundary ? 1. : phiv; // try to comment this LinearCG(newLaplace2phi,phiv,rhsv,nbiter=1000); cout << phiv.l2 << endl;