-- 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 : //    The Laplacian operator is applied on a square with a reentrant corner.
    4 : //
    5 : //    The geometry is defined by an internal angle PI < OMEGA <= 2PI.
    6 : //
    7 : //    -uxx - uyy = f in the region.
    8 : //    u = g on the boundary.
    9 : //
   10 : //    F(X,Y) = 0
   11 : //    ALPHA = PI / OMEGA
   12 : //    R = sqrt ( X^2 + Y^2 )
   13 : //    THETA = arctan ( Y / X )
   14 : //    G(X,Y) = R^ALPHA * SIN ( ALPHA * THETA)
   15 : //
   16 : //  Location:
   17 : //
   18 : //    http://people.sc.fsu.edu/~jburkardt/freefem_src/mitchell_02/mitchell_02a.edp
   19 : //
   20 : //  Licensing:
   21 : //
   22 : //    This code is distributed under the GNU LGPL license.
   23 : //
   24 : //  Modified:
   25 : //
   26 : //    24 May 2020
   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 : cout << "\n";
   45 : cout << "mitchell_02a\n";
   46 : cout << "  FreeFem++ version\n";
   47 : cout << "  Reentrant corner\n";
   48 : 
   49 : real omega = pi + 0.01;
   50 : 
   51 : cout << "  Omega = " << omega << "\n";
   52 : int p = 5;
   53 : 
   54 : border b1 ( t = 0.0, +1.0 )   { x = t;    y = 0.0; label = 1; }
   55 : real l1 = 1.0;
   56 : int n1 = ( round ) ( l1 * p + 1 );
   57 : cout << "  N1 = " << n1 << "\n";
   58 : border b2 ( t = 0.0, +1.0 )   { x = +1.0; y = t;   label = 1; }
   59 : real l2 = 1.0;
   60 : int n2 = ( round ) ( l2 * p + 1 );
   61 : cout << "  N2 = " << n2 << "\n";
   62 : border b3 ( t = 1.0, -1.0 )   { x = t;    y = 1.0; label = 1; }
   63 : real l3 = 2.0;
   64 : int n3 = ( round ) ( l3 * p + 1 );
   65 : cout << "  N3 = " << n3 << "\n";
   66 : //
   67 : //  What we do depends on OMEGA!
   68 : //  Here, we are assuming that PI <= OMEGA <= 5 PI / 4.
   69 : //
   70 : real xc = -1.0;
   71 : real yc = - tan ( omega );
   72 : real rc = - 1.0 / cos ( omega );
   73 : cout << "  XC = " << xc << "\n";
   74 : cout << "  YC = " << yc << "\n";
   75 : cout << "  RC = " << rc << "\n";
   76 : 
   77 : border b4 ( t = 1.0, yc ) { x = xc; y = t; label = 1; }
   78 : real l4 = 1.0 - yc;
   79 : int n4 = ( round ) ( l4 * p + 1 );
   80 : cout << "  N4 = " << n4 << "\n";
   81 : border b7 ( t = rc, 0.0 ) { x = t * cos ( omega ); y = t * sin ( omega ); label = 1; }
   82 : real l7 = sqrt ( xc * xc + yc * yc );
   83 : int n7 = ( round ) ( l7 * p + 1 );
   84 : cout << "  N7 = " << n7 << "\n";
   85 : //
   86 : //  Th: the mesh.
   87 : //
   88 : mesh Th = buildmesh ( b1 ( n1 ) + b2 ( n2 ) + b3 ( n3 ) + b4 ( n4 ) + b7 ( n7 ) );
   89 : //
   90 : //  Vh: the finite element space defined over Th, using P1 basis functions.
   91 : //
   92 : fespace Vh ( Th, P1 );
   93 : //
   94 : //  Define U, V, and F, piecewise continuous functions over Th.
   95 : //
   96 : Vh u;
   97 : Vh v;
   98 : //
   99 : //  Define the right hand side function F.
  100 : //
  101 : func f = 0.0;
  102 : //
  103 : //  Define the boundary condition function G.
  104 : //
  105 : real alpha = pi / omega;
  106 : func g = pow ( x * x + y * y, alpha / 2.0 ) * sin ( alpha * atan2 ( y, x ) );
  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 * v ) 
  113 :   + on ( 1, u = g );
  114 : //
  115 : //  Plot the solution.
  116 : //
  117 : plot ( u, wait = 1, fill = true, ps = "mitchell_02a_u.eps" );
  118 : //
  119 : //  Plot the mesh.
  120 : //
  121 : plot ( Th, wait = 1, ps = "mitchell_02a_mesh.eps" );
  122 : //
  123 : //  Save the mesh file.
  124 : //
  125 : savemesh ( Th, "mitchell_02a.msh" );
  126 : //
  127 : //  Terminate.
  128 : //
  129 : cout << "\n";
  130 : cout << "mitchell_02a\n";
  131 : cout << "  Normal end of execution.\n";
  132 :  sizestack + 1024 =1944  ( 920 )


mitchell_02a
  FreeFem++ version
  Reentrant corner
  Omega = 3.15159
  N1 = 6
  N2 = 6
  N3 = 11
  XC = -1
  YC = -0.0100003
  RC = 1.00005
  N4 = 6
  N7 = 6
  --  mesh:  Nb of Triangles =    163, Nb of Vertices 100
  -- Solve : 
          min -0.0199362  max 1.00634
  number of required edges : 0

mitchell_02a
  Normal end of execution.
times: compile 0.006952s, execution 0.005043s,  mpirank:0
 CodeAlloc : nb ptr  3812,  size :482656 mpirank: 0
Ok: Normal End