-- FreeFem++ v4.14 (mer. 06 mars 2024 16:59:04 CET - git v4.14-1-g2b2052ae) file : mitchell_11.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // mitchell_11.edp 2 : // 3 : // Discussion: 4 : // 5 : // -del p del(u) = 0 in [-1,+1]x[-1,+1]. 6 : // u = g on the boundary. 7 : // 8 : // p = a1*I in Q1 and Q3 9 : // = a2*I in Q2 and Q4 10 : // f = 0 11 : // g = r^gamma * mu(theta) 12 : // 13 : // mu(theta) = 14 : // Q1: cos ( gamma * ( pi/2-sigma) ) * cos ( gamma * ( theta-pi/2+rho ) ) 15 : // Q2: cos ( gamma * rho ) * cos ( gamma * ( theta-pi+sigma ) ) 16 : // Q3: cos ( gamma * sigma ) * cos ( gamma * ( theta-pi-rho ) ) 17 : // Q4: cos ( gamma * ( pi/2 - rho) ) * cos ( gamma * ( theta-3pi/2-sigma ) ) 18 : // 19 : // Parameter values: 20 : // 21 : // a1 = 161.4476387975881 22 : // a2 = 1.0 23 : // gamma = 0.1 24 : // rho = pi / 4 25 : // sigma = -14.92256510455152 26 : // 27 : // Note that the formulation in the Mitchell reference uses different 28 : // names for variables, a somewhat different formulation for mu(theta), 29 : // some typographical errors, and does not seem to produce the result 30 : // displayed in the figure. The data and formulation from the Morin 31 : // reference seems more satisfactory. 32 : // 33 : // I CANNOT get this problem to give me anything like the results 34 : // displayed by Mitchell, nor can I explain the lack of symmetry. 35 : // 36 : // Licensing: 37 : // 38 : // This code is distributed under the MIT license. 39 : // 40 : // Modified: 41 : // 42 : // 09 February 2015 43 : // 44 : // Author: 45 : // 46 : // John Burkardt 47 : // 48 : // Reference: 49 : // 50 : // Frederic Hecht, 51 : // Freefem++, 52 : // Third Edition, version 3.22 53 : // 54 : // William Mitchell, 55 : // A collection of 2D elliptic problems for testing adaptive 56 : // grid refinement algorithms, 57 : // Applied Mathematics and Computation, 58 : // Volume 220, 1 September 2013, pages 350-364. 59 : // 60 : // Pedro Morin, Ricardo Nochetto, Kunibert Siebert, 61 : // Data oscillation and convergence of adaptive FEM, 62 : // SIAM Journal on Numerical Analysis, 63 : // Volume 38, Number 2, 2000, pages 466-488. 64 : // 65 : cout << "\n"; 66 : cout << "mitchell_11():\n"; 67 : cout << " FreeFem++ version\n"; 68 : cout << " Solve the intersecting interface problem.\n"; 69 : // 70 : // Set parameters. 71 : // 72 : real a1 = 161.4476387975881; 73 : real a2 = 1.0; 74 : real gam = 0.1; 75 : real rho = pi / 4.0; 76 : real sigma = -14.92256510455152; 77 : // 78 : // Set the borders. 79 : // 80 : border bottom ( t = -1.0, 1.0 ) { x = t; y = -1.0; label = 1; } 81 : border right ( t = -1.0, 1.0 ) { x = 1.0; y = t; label = 1; } 82 : border top ( t = 1.0, -1.0 ) { x = t; y = +1.0; label = 1; } 83 : border left ( t = 1.0, -1.0 ) { x = -1.0; y = t; label = 1; } 84 : // 85 : // Define Th, the triangulation of the "left" side of the boundaries. 86 : // 87 : int n = 20; 88 : mesh Th = buildmesh ( bottom ( n ) + right ( n ) + top ( n ) + left ( n ) ); 89 : // 90 : // Define 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 exact solution function G = r^gam * mu(theta). 100 : // 101 : //func g { 102 : func real g ( real xx, real yy ) 103 : { 104 : real r; 105 : real rg; 106 : real theta; 107 : real value; 108 : 109 : r = sqrt ( xx * xx + yy * yy ); 110 : rg = pow ( r, gam ); 111 : 112 : theta = atan2 ( yy, xx ); 113 : 114 : if ( 0.0 <= xx && 0.0 <= yy ) 115 : { 116 : value = rg * cos ( gam * ( 0.5 * pi - sigma ) ) 117 : * cos ( gam * ( theta - 0.5 * pi + rho ) ); 118 : } 119 : else if ( xx < 0.0 && 0.0 <= yy ) 120 : { 121 : value = rg * cos ( gam * rho ) 122 : * cos ( gam * ( theta - pi + sigma ) ); 123 : } 124 : else if ( xx < 0.0 && yy < 0.0 ) 125 : { 126 : value = rg * cos ( gam * sigma ) 127 : * cos ( gam * ( theta - pi - rho ) ); 128 : } 129 : else if ( 0.0 <= xx && yy < 0.0 ) 130 : { 131 : value = rg * cos ( gam * ( 0.5 * pi - rho ) ) 132 : * cos ( gam * ( theta - 1.5 * pi - sigma ) ); 133 : } 134 : else 135 : { 136 : cout << "WHAT THE HELL? X = " << xx << " Y = " << yy << "\n"; 137 : } 138 : return value; 139 : } 140 : // 141 : // Define right hand side function F. 142 : // 143 : func f = 0.0; 144 : // 145 : // Define coefficient function P. 146 : // 147 : func real p ( real xx, real yy ) 148 : { 149 : real value; 150 : if ( 0.0 <= xx * yy ) 151 : { 152 : value = a1; 153 : } 154 : else 155 : { 156 : value = a2; 157 : } 158 : return value; 159 : } 160 : // 161 : // Solve the variational problem. 162 : // 163 : solve Laplace ( u, v ) 164 : = int2d ( Th ) ( p(x,y) * ( dx(u)*dx(v) + dy(u)*dy(v) ) ) 165 : - int2d ( Th ) ( f * v ) 166 : + on ( 1, u = g(x,y) ); 167 : // 168 : // Plot the solution. 169 : // 170 : plot ( u, wait = true, fill = true, ps = "mitchell_11_u.ps" ); 171 : // 172 : // Plot the mesh. 173 : // 174 : plot ( Th, wait = true, ps = "mitchell_11_mesh.ps" ); 175 : // 176 : // Save the mesh file. 177 : // 178 : savemesh ( Th, "mitchell_11.msh" ); 179 : // 180 : // Terminate. 181 : // 182 : cout << "\n"; 183 : cout << "mitchell_11():\n"; 184 : cout << " Normal end of execution.\n"; 185 : 186 : exit ( 0 ); 187 : sizestack + 1024 =1904 ( 880 ) mitchell_11(): FreeFem++ version Solve the intersecting interface problem. -- mesh: Nb of Triangles = 938, Nb of Vertices 510 -- Solve : min -0.0812259 max 0.647446 number of required edges : 0 mitchell_11(): Normal end of execution. current line = 186 exit(0) err code 0 , mpirank 0 CodeAlloc : nb ptr 4159, size :530696 mpirank: 0 Ok: Normal End