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