-- FreeFem++ v  3.320001 (date Ven  7 nov 2014 14:28:55 CET)
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //  Discussion:
    2 : //
    3 : //    Oscillatory
    4 : //    -Laplace(u) -u/(alpha+r)^4 = f on the unit square.
    5 : //    r = sqrt ( x^2 + y^2 )
    6 : //    u = u0 on the boundary.
    7 : //    u0 = sin(1/(alpha+r))
    8 : //    f = -Laplace(u0) -u0/(alpha+r)^4
    9 : //
   10 : //    The parameter alpha affects the number of oscillations in the solution.
   11 : //
   12 : //    Suggested parameter values:
   13 : //    * alpha = 1/(10*pi) is relatively easy.
   14 : //    * alpha = 1/(50*pi) is harder.
   15 : //
   16 : //  Location:
   17 : //
   18 : //    http://people.sc.fsu.edu/~jburkardt/examples/mitchell_freefem++/test08b.edp
   19 : //
   20 : //  Licensing:
   21 : //
   22 : //    This code is distributed under the GNU LGPL license.
   23 : //
   24 : //  Modified:
   25 : //
   26 : //    20 December 2014
   27 : //
   28 : //  Author:
   29 : //
   30 : //    John Burkardt
   31 : //
   32 : //  Reference:
   33 : //
   34 : //    Frederic Hecht,
   35 : //    Freefem++,
   36 : //    Third Edition, version 3.22
   37 : //
   38 : //    William Mitchell,
   39 : //    A collection of 2D elliptic problems for testing adaptive
   40 : //    grid refinement algorithms,
   41 : //    Applied Mathematics and Computation,
   42 : //    Volume 220, 1 September 2013, pages 350-364.
   43 : //
   44 : border bottom ( t = 0.0, 1.0 )   { x = t;   y = 0.0;  label = 1; }
   45 : border right  ( t = 0.0, 1.0 )   { x = 1.0; y = t;    label = 1; }
   46 : border top    ( t = 1.0, 0.0 )   { x = t;   y = 1.0;  label = 1; }
   47 : border left   ( t = 1.0, 0.0 )   { x = 0.0; y = t;    label = 1; }
   48 : //
   49 : //  Define Th, the triangulation of the "left" side of the boundaries.
   50 : //
   51 : int n = 30;
   52 : mesh Th = buildmesh ( bottom ( n ) + right ( n ) + top ( n ) + left ( n ) );
   53 : //
   54 : //  Define Vh, the finite element space defined over Th, using P1 basis functions.
   55 : //
   56 : fespace Vh ( Th, P1 );
   57 : //
   58 : //  Define U, V, and F, piecewise continuous functions over Th.
   59 : //
   60 : Vh u;
   61 : Vh v;
   62 : //
   63 : //  Define problem parameters.
   64 : //
   65 : real alpha = 1.0 / 50.0 / pi;
   66 : //
   67 : //  Define the right hand side function F.
   68 : //
   69 : func f = 
   70 :   ( 
   71 :     ( - alpha * alpha + x * x + y * y ) 
   72 :     * cos ( 1.0 / ( alpha + sqrt ( x * x + y * y ) ) ) 
   73 :     - ( x * x + y * y ) 
   74 :     * sin ( 1.0 / ( alpha + sqrt ( x * x + y * y ) ) )
   75 :   ) 
   76 :   / sqrt ( x * x + y * y ) 
   77 :   / pow ( alpha + sqrt ( x * x + y * y ), 4 );
   78 : //
   79 : //  Define the boundary function G.
   80 : //
   81 : func g = sin ( 1.0 / alpha + sqrt ( x * x + y * y ) );
   82 : //
   83 : //  Solve the variational problem.
   84 : //
   85 : solve Laplace ( u, v ) 
   86 :   = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
   87 :   - int2d ( Th ) ( u * v / pow ( alpha + sqrt ( x * x + y * y ), 4 ) )
   88 :   - int2d ( Th ) ( f * v ) 
   89 :   + on ( 1, u = g );
   90 : //
   91 : //  Plot the solution.
   92 : //
   93 : plot ( u, wait = 1, fill = true, ps = "test08b_u.eps" );
   94 : //
   95 : //  Plot the mesh.
   96 : //
   97 : plot ( Th, wait = 1, ps = "test08b_mesh.eps" );
   98 : //
   99 : //  Save the mesh file.
  100 : //
  101 : savemesh ( Th, "test08b.msh" );
  102 : 
  103 :  sizestack + 1024 =5336  ( 4312 )

  --  mesh:  Nb of Triangles =   2096, Nb of Vertices 1109
  -- Solve : 
          min -1.39905  max 0.987766
  number of required edges : 0
times: compile 0.01s, execution 0.08s,  mpirank:0
 CodeAlloc : nb ptr  2768,  size :346520 mpirank: 0
Ok: Normal End