-- 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++/test02e.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 : 47 : // 48 : // We can't use omega = 2pi, or even 2pi - 0.01. 49 : // 50 : real omega = 8.0 * pi / 4.0 - 0.05; 51 : 52 : cout << " Omega = " << omega << "\n"; 53 : int p = 5; 54 : 55 : border b1 ( t = 0.0, +1.0 ) { x = t; y = 0.0; label = 1; } 56 : real l1 = 1.0; 57 : int n1 = ( round ) ( l1 * p + 1 ); 58 : cout << " N1 = " << n1 << "\n"; 59 : border b2 ( t = 0.0, +1.0 ) { x = +1.0; y = t; label = 1; } 60 : real l2 = 1.0; 61 : int n2 = ( round ) ( l2 * p + 1 ); 62 : cout << " N2 = " << n2 << "\n"; 63 : border b3 ( t = 1.0, -1.0 ) { x = t; y = 1.0; label = 1; } 64 : real l3 = 2.0; 65 : int n3 = ( round ) ( l3 * p + 1 ); 66 : cout << " N3 = " << n3 << "\n"; 67 : border b4 ( t = 1.0, -1.0 ) { x = -1.0; y = t; label = 1; } 68 : real l4 = 2.0; 69 : int n4 = ( round ) ( l4 * p + 1 ); 70 : cout << " N4 = " << n4 << "\n"; 71 : border b5 ( t = -1.0, +1.0 ) { x = t; y = -1.0; label = 1; } 72 : real l5 = 2.0; 73 : int n5 = ( round ) ( l5 * p + 1 ); 74 : cout << " N5 = " << n5 << "\n"; 75 : real xc = 1.0; 76 : real yc = tan ( omega ); 77 : real rc = 1.0 / cos ( omega ); 78 : 79 : border b6 ( t = -1.0, yc ) { x = 1.0; y = t; label = 1; } 80 : real l6 = yc + 1.0; 81 : int n6 = ( round ) ( l6 * p + 1 ); 82 : cout << " N6 = " << n6 << "\n"; 83 : border b7 ( t = rc, 0.0 ) { x = t * cos ( omega ); y = t * sin ( omega ); label = 1; } 84 : real l7 = rc; 85 : int n7 = ( round ) ( l7 * p + 1 ); 86 : cout << " N7 = " << n7 << "\n"; 87 : // 88 : // Create the mesh. 89 : // 90 : mesh Th = buildmesh ( b1 ( n1 ) + b2 ( n2 ) + b3 ( n3 ) + b4 ( n4 ) 91 : + b5 ( n5 ) + b6 ( n6 ) + b7 ( n7 ) ); 92 : // 93 : // Define Vh, the finite element space defined over Th, using P1 basis functions. 94 : // 95 : fespace Vh ( Th, P1 ); 96 : // 97 : // Define U, V, and F, piecewise continuous functions over Th. 98 : // 99 : Vh u; 100 : Vh v; 101 : // 102 : // Define the right hand side function F. 103 : // 104 : func f = 0.0; 105 : // 106 : // Define the boundary condition function G. 107 : // 108 : real alpha = pi / omega; 109 : func g = pow ( x * x + y * y, alpha / 2.0 ) * sin ( alpha * atan2 ( y, x ) ); 110 : // 111 : // Solve the variational problem. 112 : // 113 : solve Laplace ( u, v ) 114 : = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) ) 115 : - int2d ( Th ) ( f * v ) 116 : + on ( 1, u = g ); 117 : // 118 : // Plot the solution. 119 : // 120 : plot ( u, wait = 1, fill = true, ps = "test02e_u.eps" ); 121 : // 122 : // Plot the mesh. 123 : // 124 : plot ( Th, wait = 1, ps = "test02e_mesh.eps" ); 125 : // 126 : // Save the mesh file. 127 : // 128 : savemesh ( Th, "test02e.msh" ); 129 : 130 : sizestack + 1024 =2008 ( 984 ) Omega = 6.23319 N1 = 6 N2 = 6 N3 = 11 N4 = 11 N5 = 11 N6 = 6 N7 = 6 -- mesh: Nb of Triangles = 291, Nb of Vertices 174 -- Solve : min -1.10447 max 1.10447 number of required edges : 0 Some Double edge in the mesh, the number is 57 nbe4=56 Fatal error in the meshgenerator 1002 current line = 128 Meshing error: Bamg number : 1002, Meshing error: Bamg number : 1002, err code 4 , mpirank 0