-- 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 : //    -uxx-uyy = f on [-1,+1]x[-1,+1]
    4 : //    u = g on the boundary.
    5 : //
    6 : //    f = -gxx-gyy.
    7 : //    g = cos(pi*y/2)                        for x <= beta*(y+1)
    8 : //      = cos(pi*y/2)+(x-beta*(y+1))**alpha  for      beta*(y+1) < x
    9 : //
   10 : //    Suggested parameter values:
   11 : //    * alpha = 2.5, beta = 0.0
   12 : //    * alpha = 1.1, beta = 0.0
   13 : //    * alpha = 1.5, beta = 0.6
   14 : //
   15 : //  Location:
   16 : //
   17 : //    http://people.sc.fsu.edu/~jburkardt/examples/mitchell_freefem++/test10a.edp
   18 : //
   19 : //  Licensing:
   20 : //
   21 : //    This code is distributed under the GNU LGPL license.
   22 : //
   23 : //  Modified:
   24 : //
   25 : //    23 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 : 
   44 : //
   45 : //  Set parameters.
   46 : //
   47 : real alpha = 2.5;
   48 : real beta = 0.0;
   49 : //
   50 : //  Set the borders.
   51 : //
   52 : border bottom ( t = -1.0,  1.0 )  { x = t;     y = -1.0; label = 1; }
   53 : border right  ( t = -1.0,  1.0 )  { x = 1.0;   y = t;    label = 1; }
   54 : border top    ( t =  1.0, -1.0 )  { x = t;     y = +1.0; label = 1; }
   55 : border left   ( t =  1.0, -1.0 )  { x = -1.0;  y = t;    label = 1; }
   56 : //
   57 : //  Define Th, the triangulation of the "left" side of the boundaries.
   58 : //
   59 : int n = 10;
   60 : mesh Th = buildmesh ( bottom ( n ) + right ( n ) + top ( n ) + left ( n ) );
   61 : //
   62 : //  Define Vh, the finite element space defined over Th, using P1 basis functions.
   63 : //
   64 : fespace Vh ( Th, P1 );
   65 : //
   66 : //  Define U, V, and F, piecewise continuous functions over Th.
   67 : //
   68 : Vh u;
   69 : Vh v;
   70 : //
   71 : //  Define function G.
   72 : //
   73 : func real g ( real alpha, real beta )
   74 : {
   75 :   real value;
   76 :   if ( x <= beta * ( y + 1.0 ) )
   77 :   {
   78 :     value = cos ( pi * y / 2.0 );
   79 :   }
   80 :   else
   81 :   {
   82 :     value = cos ( pi * y / 2.0 ) + pow ( x - beta * ( y + 1.0 ), alpha );
   83 :   }
   84 :   return value;
   85 : }
   86 : //
   87 : //  Define function F.
   88 : //
   89 : func real f ( real alpha, real beta )
   90 : {
   91 :   real value;
   92 :   if ( x <= beta * ( y + 1.0 ) )
   93 :   {
   94 :     value = 0.25 * pi * pi * cos ( pi * y / 2.0 );
   95 :   }
   96 :   else
   97 :   {
   98 :     value = 0.25 * pi * pi * cos ( pi * y / 2.0 ) 
   99 :       + alpha * ( alpha - 1 ) * ( 1.0 + beta * beta )
  100 :       * pow ( x - beta * ( y + 1.0 ), alpha - 2 );
  101 :   }
  102 :   return value;
  103 : }
  104 : //
  105 : //  Solve the variational problem.
  106 : //
  107 : solve Laplace ( u, v ) 
  108 :   = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
  109 :   - int2d ( Th ) ( f ( alpha, beta ) * v ) 
  110 :   + on ( 1, u = g ( alpha, beta ) );
  111 : //
  112 : //  Plot the solution.
  113 : //
  114 : plot ( u, wait = 1, fill = true, ps = "test10a_u.eps" );
  115 : //
  116 : //  Plot the mesh.
  117 : //
  118 : plot ( Th, wait = 1, ps = "test10a_mesh.eps" );
  119 : //
  120 : //  Save the mesh file.
  121 : //
  122 : savemesh ( Th, "test10a.msh" );
  123 : 
  124 :  sizestack + 1024 =1856  ( 832 )

  --  mesh:  Nb of Triangles =    240, Nb of Vertices 141
  -- Solve : 
          min 6.12323e-17  max 2.04549
  number of required edges : 0
times: compile 0s, execution 0.04s,  mpirank:0
 CodeAlloc : nb ptr  2740,  size :345072 mpirank: 0
Ok: Normal End