// -*- FreeFem++ -*- // Time-stamp: "2020-07-18 21:40:41 fujiwara" // Original : https://modules.freefem.org/modules/poisson/ // Parameters func f = 1.; // Mesh int nn = 2; //Mesh quality mesh Th = square(nn, nn); // Fespace //func Pk = P1nc; // equivalent load "Element_CR" func Pk = CR2d; fespace Uh(Th, Pk); Uh u; // Macro macro grad(A) [dx(A), dy(A)] // // Problem varf vPoisson (u, uh) = int2d(Th)( grad(u)' * grad(uh) ) + int2d(Th)( f * uh ) + on(1, 2, 3, 4, u=0) ; matrix Poisson = vPoisson(Uh, Uh, solver=sparsesolver); real[int] PoissonBoundary = vPoisson(0, Uh); u[] = Poisson^-1 * PoissonBoundary; // Plot plot(u, nbiso=30, fill=true, value=true, cmm="A"); //---------------------------------------------------------------------- // output values at each middle point on edges //---------------------------------------------------------------------- ofstream fout("Poisson2d-u-CR.dat"); for(int it = 0; it < Th.nt; it++){ real v0x = Th( Th[it][0] ).x; real v0y = Th( Th[it][0] ).y; real v1x = Th( Th[it][1] ).x; real v1y = Th( Th[it][1] ).y; real v2x = Th( Th[it][2] ).x; real v2y = Th( Th[it][2] ).y; // middle points on each edge real e0x = ( v1x + v2x ) / 2; real e0y = ( v1y + v2y ) / 2; real e1x = ( v2x + v0x ) / 2; real e1y = ( v2y + v0y ) / 2; real e2x = ( v0x + v1x ) / 2; real e2y = ( v0y + v1y ) / 2; real u0 = u(e0x,e0y); real u1 = u(e1x,e1y); real u2 = u(e2x,e2y); fout << "# Triangle " << it << endl; fout << e0x << " " << e0y << " " << u0 << endl; fout << e1x << " " << e1y << " " << u1 << endl; fout << e2x << " " << e2y << " " << u2 << endl; fout << endl; } //---------------------------------------------------------------------- // End of file //----------------------------------------------------------------------