-- FreeFem++ v4.6 (Thu Apr  2 15:47:38 CEST 2020 - git v4.6)
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //  contour.edp
    2 : //
    3 : //  Discussion:
    4 : //
    5 : //    Let an ellipse have semimajor axis 2 and semiminor axis 1:
    6 : //      (x/2)^2 + (y/1)^2 = 1.
    7 : //    Let Gamma1 be the ellipse boundary from 0 to 4pi/3, and
    8 : //    Gamm2 the boundary from 4pi/3 to 2pi.
    9 : //
   10 : //    Solve:
   11 : //      -uxx - uyy = f in Gamma1 + Gamma2,
   12 : //      u = x on Gamma1.
   13 : //      du/dn = 0, natural boundary condition on Gamma2.
   14 : //
   15 : //  Location:
   16 : //
   17 : //    http://people.sc.fsu.edu/~jburkardt/freefem_src/contour/contour.edp
   18 : //
   19 : //  Licensing:
   20 : //
   21 : //    This code is distributed under the GNU LGPL license.
   22 : //
   23 : //  Modified:
   24 : //
   25 : //    23 May 2020
   26 : //
   27 : //  Reference:
   28 : //
   29 : //    Frederic Hecht,
   30 : //    Freefem++,
   31 : //    Third Edition, version 3.22
   32 : //
   33 : cout << "\n";
   34 : cout << "gnuplot_contour:\n";
   35 : cout << "  FreeFem++ version\n";
   36 : cout << "  Demonstrate how to write data that GNUPLOT can 
  ... : use\n";
   37 : cout << "  to make a contour plot of the solution.\n";
   38 : //
   39 : //  Define Gamma1 and Gamma2, the boundaries.
   40 : //  A and B are the lengths of the semimajor and semiminor axes:
   41 : //  THETA specifies the angle at which we switch from GAMMA1 to GAMMA2.
   42 : //
   43 : real a = 2.0;
   44 : real b = 1.0;
   45 : real theta = 4.0 * pi / 3.0;
   46 : border Gamma1 ( t = 0, theta )    { x = a * cos(t); y = b*sin(t); }
   47 : border Gamma2 ( t = theta, 2*pi ) { x = a * cos(t); y = b*sin(t); }
   48 : //
   49 : //  Th: the triangulation of the "left" side of the boundaries.
   50 : //
   51 : mesh Th = buildmesh ( Gamma1(100) + Gamma2(50) );
   52 : //
   53 : //  Plot the mesh.
   54 : //
   55 : plot ( Th, wait = 0, ps = "gnuplot_contour_mesh.ps" );
   56 : //
   57 : //  Define Vh, the finite element space defined over Th, using P2 basis functions.
   58 : //
   59 : fespace Vh ( Th, P2 );
   60 : //
   61 : //  Define U, V, and F, piecewise P2 continuous functions over Th.
   62 : //  F is a constant function.
   63 : //
   64 : Vh u;
   65 : Vh v;
   66 : Vh f = 1.0;
   67 : //
   68 : //  Define G, a function for the Dirichlet boundary conditions.
   69 : //
   70 : func g = x;
   71 : //
   72 : //  Define the weak form of the PDE.
   73 : //
   74 : problem poisson ( u, v ) 
   75 :   = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
   76 :   - int2d ( Th ) ( f * v ) 
   77 :   + on ( Gamma1, u = g );
   78 : //
   79 : //  Solve the PDE.
   80 : //
   81 : poisson;
   82 : //
   83 : //  Write gnuplot data.
   84 : //
   85 : ofstream dataunit ( "gnuplot_contour_data.txt" );
   86 : for ( int i = 0; i < Th.nt; i++ )
   87 : {
   88 :   for ( int j = 0; j <= 3; j++ )
   89 :   {
   90 :     int jj = ( j % 3 );
   91 :     dataunit << "  " << Th[i][jj].x 
   92 :              << "  " << Th[i][jj].y
   93 :              << "  " << u[][Vh(i,jj)] << "\n";
   94 :   }
   95 : }
   96 : //
   97 : //  Write gnuplot commands.
   98 : //
   99 : ofstream commandunit ( "gnuplot_contour_commands.txt" );
  100 : commandunit << "# gnuplot_contour_commands.txt\n";
  101 : commandunit << "# usage: gnuplot < gnuplot_contour_commands.txt\n
  ... : ";
  102 : commandunit << "#\n";
  103 : commandunit << "set term png\n";
  104 : commandunit << "set size ratio -1\n";
  105 : commandunit << "set output 'gnuplot_contour.png'\n";
  106 : commandunit << "set pm3d\n";
  107 : commandunit << "unset surface\n";
  108 : commandunit << "set view map\n";
  109 : commandunit << "set contour\n";
  110 : commandunit << "set nokey\n";
  111 : commandunit << "set cntrparam cubicspline  # smooth out the lines
  ... : \n";
  112 : commandunit << "set cntrparam levels 50    # sets the num of cont
  ... : our lines\n";
  113 : commandunit << "set pm3d interpolate 20,20 # interpolate the colo
  ... : r\n";
  114 : commandunit << "set palette model RGB defined ( 0'black', 1'blue'
  ... : , 2'cyan',3'green',4'yellow',5'red',8'purple' )\n";
  115 : commandunit << "set xlabel '<--X-->'\n";
  116 : commandunit << "set ylabel '<--Y-->'\n";
  117 : commandunit << "set title 'Solution contours'\n";
  118 : commandunit << "set timestamp\n";
  119 : commandunit << "set dgrid3d 10, 10, 1\n";
  120 : commandunit << "splot 'gnuplot_contour_data.txt' with lines lt 1\n";
  121 : //
  122 : //  Terminate.
  123 : //
  124 : cout << "\n";
  125 : cout << "gnuplot_contour:\n";
  126 : cout << "  Normal end of execution.\n";
  127 : 
  128 : exit ( 0 );
  129 :  sizestack + 1024 =2464  ( 1440 )


gnuplot_contour:
  FreeFem++ version
  Demonstrate how to write data that GNUPLOT can use
  to make a contour plot of the solution.
  --  mesh:  Nb of Triangles =   2568, Nb of Vertices 1360
  -- Solve : 
          min -2  max 2

gnuplot_contour:
  Normal end of execution.
  current line = 128
exit(0)
 err code 0 ,  mpirank 0
 CodeAlloc : nb ptr  3730,  size :479032 mpirank: 0
Ok: Normal End