-- FreeFem++ v4.6 (Thu Apr  2 15:47:38 CEST 2020 - git v4.6)
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //  Discussion:
    2 : //
    3 : //    - uxx - uyy = f on the unit square.
    4 : //    u = g on the boundary.
    5 : //
    6 : //    f = - gxx - gyy
    7 : //    g = atan ( alpha * (sqrt((x-xc)^2+(y-yc)^2)-r0) )
    8 : //
    9 : //    Suggested parameter values:
   10 : //    * alpha =   20, xc = -0.05, yc = -0.05, r0 = 0.7
   11 : //    * alpha = 1000, xc = -0.05, yc = -0.05, r0 = 0.7
   12 : //    * alpha = 1000, xc =  1.5,  yc =  0.25, r0 = 0.92
   13 : //    * alpha =   50, xc =  0.5,  yc =  0.5,  r0 = 0.25
   14 : //
   15 : //  Location:
   16 : //
   17 : //    http://people.sc.fsu.edu/~jburkardt/freefem_src/mitchell_09/mitchell_09a.edp
   18 : //
   19 : //  Licensing:
   20 : //
   21 : //    This code is distributed under the GNU LGPL license.
   22 : //
   23 : //  Modified:
   24 : //
   25 : //    22 December 2014
   26 : //
   27 : //  Author:
   28 : //
   29 : //    John Burkardt
   30 : //
   31 : //  Reference:
   32 : //
   33 : //    Frederic Hecht,
   34 : //    Freefem++,
   35 : //    Third Edition, version 3.22
   36 : //
   37 : //    William Mitchell,
   38 : //    A collection of 2D elliptic problems for testing adaptive
   39 : //    grid refinement algorithms,
   40 : //    Applied Mathematics and Computation,
   41 : //    Volume 220, 1 September 2013, pages 350-364.
   42 : //
   43 : cout << "\n";
   44 : cout << "mitchell_09a\n";
   45 : cout << "  FreeFem++ version\n";
   46 : cout << "  The wave front problem.\n";
   47 : 
   48 : border bottom ( t = 0.0, 1.0 )   { x = t;   y = 0.0;  label = 1; }
   49 : border right  ( t = 0.0, 1.0 )   { x = 1.0; y = t;    label = 1; }
   50 : border top    ( t = 1.0, 0.0 )   { x = t;   y = 1.0;  label = 1; }
   51 : border left   ( t = 1.0, 0.0 )   { x = 0.0; y = t;    label = 1; }
   52 : //
   53 : //  Define Th, the triangulation of the "left" side of the boundaries.
   54 : //
   55 : int n = 10;
   56 : mesh Th = buildmesh ( bottom ( n ) + right ( n ) + top ( n ) + left ( n ) );
   57 : //
   58 : //  Define Vh, the finite element space defined over Th, using P1 basis functions.
   59 : //
   60 : fespace Vh ( Th, P1 );
   61 : //
   62 : //  Define U, V, and F, piecewise continuous functions over Th.
   63 : //
   64 : Vh u;
   65 : Vh v;
   66 : //
   67 : //  Define problem parameters.
   68 : //
   69 : real alpha = 20.0;
   70 : real xc = - 0.05;
   71 : real yc = - 0.05;
   72 : real r0 = 0.7;
   73 : //
   74 : //  Define the right hand side function F.
   75 : //
   76 : func f = 
   77 :   ( 
   78 :     - alpha + pow ( alpha, 3 ) 
   79 :     * ( - r0 * r0 + pow ( x - xc, 2 ) + pow ( y - yc, 2  ) )
   80 :   )
   81 :   /
   82 :   ( 
   83 :     sqrt ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) 
   84 :     * pow ( 
   85 :       -1.0 + alpha * alpha * 
   86 :       ( 
   87 :         - r0 * r0 - pow ( x - xc, 2 ) - pow ( y - yc, 2 )
   88 :         + sqrt ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) 
   89 :       )
   90 :     , 2 ) 
   91 :   );
   92 : //
   93 : //  Define the boundary function G.
   94 : //
   95 : func g = atan ( alpha * ( sqrt ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) - r0 ) );
   96 : //
   97 : //  Solve the variational problem.
   98 : //
   99 : solve Laplace ( u, v ) 
  100 :   = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
  101 :   - int2d ( Th ) ( f * v ) 
  102 :   + on ( 1, u = g );
  103 : //
  104 : //  Plot the solution.
  105 : //
  106 : plot ( u, wait = 1, fill = true, ps = "mitchell_09a_u.ps" );
  107 : //
  108 : //  Plot the mesh.
  109 : //
  110 : plot ( Th, wait = 1, ps = "mitchell_09a_mesh.ps" );
  111 : //
  112 : //  Save the mesh file.
  113 : //
  114 : savemesh ( Th, "mitchell_09a.msh" );
  115 : //
  116 : //  Terminate.
  117 : //
  118 : cout << "\n";
  119 : cout << "mitchell_09a\n";
  120 : cout << "  Normal end of execution.\n";
  121 :  sizestack + 1024 =4800  ( 3776 )


mitchell_09a
  FreeFem++ version
  The wave front problem.
  --  mesh:  Nb of Triangles =    240, Nb of Vertices 141
  -- Solve : 
          min -1.49151  max 1.50718
  number of required edges : 0

mitchell_09a
  Normal end of execution.
times: compile 0.011701s, execution 0.008888s,  mpirank:0
 CodeAlloc : nb ptr  3737,  size :480880 mpirank: 0
Ok: Normal End