-- 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 [-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/freefem_src/mitchell_10/mitchell_10c.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 : cout << "\n";
   44 : cout << "mitchell_10c\n";
   45 : cout << "  FreeFem++ version\n";
   46 : cout << "  Solve the interior line singularity problem.\n";
   47 : //
   48 : //  Set parameters.
   49 : //
   50 : real alpha = 1.5;
   51 : real beta = 0.6;
   52 : //
   53 : //  Set the borders.
   54 : //
   55 : border bottom ( t = -1.0,  1.0 )  { x = t;     y = -1.0; label = 1; }
   56 : border right  ( t = -1.0,  1.0 )  { x = 1.0;   y = t;    label = 1; }
   57 : border top    ( t =  1.0, -1.0 )  { x = t;     y = +1.0; label = 1; }
   58 : border left   ( t =  1.0, -1.0 )  { x = -1.0;  y = t;    label = 1; }
   59 : //
   60 : //  Define Th, the triangulation of the "left" side of the boundaries.
   61 : //
   62 : int n = 10;
   63 : mesh Th = buildmesh ( bottom ( n ) + right ( n ) + top ( n ) + left ( n ) );
   64 : //
   65 : //  Define Vh, the finite element space defined over Th, using P1 basis functions.
   66 : //
   67 : fespace Vh ( Th, P1 );
   68 : //
   69 : //  Define U, V, and F, piecewise continuous functions over Th.
   70 : //
   71 : Vh u;
   72 : Vh v;
   73 : //
   74 : //  Define function G.
   75 : //
   76 : func real g ( real alpha, real beta )
   77 : {
   78 :   real value;
   79 :   if ( x <= beta * ( y + 1.0 ) )
   80 :   {
   81 :     value = cos ( pi * y / 2.0 );
   82 :   }
   83 :   else
   84 :   {
   85 :     value = cos ( pi * y / 2.0 ) + pow ( x - beta * ( y + 1.0 ), alpha );
   86 :   }
   87 :   return value;
   88 : }
   89 : //
   90 : //  Define function F.
   91 : //
   92 : func real f ( real alpha, real beta )
   93 : {
   94 :   real value;
   95 :   if ( x <= beta * ( y + 1.0 ) )
   96 :   {
   97 :     value = 0.25 * pi * pi * cos ( pi * y / 2.0 );
   98 :   }
   99 :   else
  100 :   {
  101 :     value = 0.25 * pi * pi * cos ( pi * y / 2.0 ) 
  102 :       + alpha * ( alpha - 1 ) * ( 1.0 + beta * beta )
  103 :       * pow ( x - beta * ( y + 1.0 ), alpha - 2 );
  104 :   }
  105 :   return value;
  106 : }
  107 : //
  108 : //  Solve the variational problem.
  109 : //
  110 : solve Laplace ( u, v ) 
  111 :   = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
  112 :   - int2d ( Th ) ( f ( alpha, beta ) * v ) 
  113 :   + on ( 1, u = g ( alpha, beta ) );
  114 : //
  115 : //  Plot the solution.
  116 : //
  117 : plot ( u, wait = 1, fill = true, ps = "mitchell_10c_u.ps" );
  118 : //
  119 : //  Plot the mesh.
  120 : //
  121 : plot ( Th, wait = 1, ps = "mitchell_10c_mesh.ps" );
  122 : //
  123 : //  Save the mesh file.
  124 : //
  125 : savemesh ( Th, "mitchell_10c.msh" );
  126 : //
  127 : //  Terminate.
  128 : //
  129 : cout << "\n";
  130 : cout << "mitchell_10c\n";
  131 : cout << "  Normal end of execution.\n";
  132 :  sizestack + 1024 =1856  ( 832 )


mitchell_10c
  FreeFem++ version
  Solve the interior line singularity problem.
  --  mesh:  Nb of Triangles =    240, Nb of Vertices 141
  -- Solve : 
          min 6.12323e-17  max 1.42332
  number of required edges : 0

mitchell_10c
  Normal end of execution.
times: compile 0.009748s, execution 0.032442s,  mpirank:0
 CodeAlloc : nb ptr  3715,  size :480712 mpirank: 0
Ok: Normal End