-- 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 : // 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 : // ALPHA = PI / OMEGA 8 : // R = sqrt ( X^2 + Y^2 ) 9 : // THETA = arctan ( Y / X ) 10 : // U(X,Y) = R^ALPHA * SIN ( ALPHA * THETA) 11 : // 12 : // Dirichlet boundary conditions. 13 : // 14 : // Surprisingly, U(X,Y) satisfies the equation: 15 : // 16 : // - Uxx - Uyy = 0 17 : // 18 : // Location: 19 : // 20 : // http://people.sc.fsu.edu/~jburkardt/examples/mitchell_freefem++/test02a.edp 21 : // 22 : // Licensing: 23 : // 24 : // This code is distributed under the GNU LGPL license. 25 : // 26 : // Modified: 27 : // 28 : // 17 December 2014 29 : // 30 : // Author: 31 : // 32 : // John Burkardt 33 : // 34 : // Reference: 35 : // 36 : // Frederic Hecht, 37 : // Freefem++, 38 : // Third Edition, version 3.22 39 : // 40 : // William Mitchell, 41 : // A collection of 2D elliptic problems for testing adaptive 42 : // grid refinement algorithms, 43 : // Applied Mathematics and Computation, 44 : // Volume 220, 1 September 2013, pages 350-364. 45 : // 46 : real omega = pi + 0.01; 47 : 48 : cout << " Omega = " << omega << "\n"; 49 : int p = 5; 50 : 51 : border b1 ( t = 0.0, +1.0 ) { x = t; y = 0.0; label = 1; } 52 : real l1 = 1.0; 53 : int n1 = ( round ) ( l1 * p + 1 ); 54 : cout << " N1 = " << n1 << "\n"; 55 : border b2 ( t = 0.0, +1.0 ) { x = +1.0; y = t; label = 1; } 56 : real l2 = 1.0; 57 : int n2 = ( round ) ( l2 * p + 1 ); 58 : cout << " N2 = " << n2 << "\n"; 59 : border b3 ( t = 1.0, -1.0 ) { x = t; y = 1.0; label = 1; } 60 : real l3 = 2.0; 61 : int n3 = ( round ) ( l3 * p + 1 ); 62 : cout << " N3 = " << n3 << "\n"; 63 : // 64 : // What we do depends on OMEGA! 65 : // Here, we are assuming that PI <= OMEGA <= 5 PI / 4. 66 : // 67 : real xc = -1.0; 68 : real yc = - tan ( omega ); 69 : real rc = - 1.0 / cos ( omega ); 70 : cout << " XC = " << xc << "\n"; 71 : cout << " YC = " << yc << "\n"; 72 : cout << " RC = " << rc << "\n"; 73 : 74 : border b4 ( t = 1.0, yc ) { x = xc; y = t; label = 1; } 75 : real l4 = 1.0 - yc; 76 : int n4 = ( round ) ( l4 * p + 1 ); 77 : cout << " N4 = " << n4 << "\n"; 78 : border b7 ( t = rc, 0.0 ) { x = t * cos ( omega ); y = t * sin ( omega ); label = 1; } 79 : real l7 = sqrt ( xc * xc + yc * yc ); 80 : int n7 = ( round ) ( l7 * p + 1 ); 81 : cout << " N7 = " << n7 << "\n"; 82 : // 83 : // Build the mesh. 84 : // 85 : mesh Th = buildmesh ( b1 ( n1 ) + b2 ( n2 ) + b3 ( n3 ) + b4 ( n4 ) + b7 ( n7 ) ); 86 : // 87 : // Define Vh, the finite element space defined over Th, using P1 basis functions. 88 : // 89 : fespace Vh ( Th, P1 ); 90 : // 91 : // Define U, V, and F, piecewise continuous functions over Th. 92 : // 93 : Vh u; 94 : Vh v; 95 : // 96 : // Define the right hand side function F. 97 : // 98 : func f = 0.0; 99 : // 100 : // Define the boundary condition function G. 101 : // 102 : real alpha = pi / omega; 103 : func g = pow ( x * x + y * y, alpha / 2.0 ) * sin ( alpha * atan2 ( y, x ) ); 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 * v ) 110 : + on ( 1, u = g ); 111 : // 112 : // Plot the solution. 113 : // 114 : plot ( u, wait = 1, fill = true, ps = "test02a_u.eps" ); 115 : // 116 : // Plot the mesh. 117 : // 118 : plot ( Th, wait = 1, ps = "test02a_mesh.eps" ); 119 : // 120 : // Save the mesh file. 121 : // 122 : savemesh ( Th, "test02a.msh" ); 123 : 124 : sizestack + 1024 =1944 ( 920 ) 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 times: compile 0s, execution 0.01s, mpirank:0 CodeAlloc : nb ptr 2833, size :346888 mpirank: 0 Ok: Normal End