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