-- 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 : // Peak 4 : // -Laplace(u) = f on the unit square. 5 : // u = u0 on the boundary. 6 : // u0 = exp(-alpha*(x-xc)^2+(y-yc)^2) 7 : // f = -(u0'') 8 : // 9 : // Discussion: 10 : // 11 : // Suggested parameter values: 12 : // * alpha = 1000, (xc,yc) = (0.5,0.5) 13 : // * alpha = 100000, (xc,yc) = (0.51,0.117) 14 : // 15 : // Location: 16 : // 17 : // http://people.sc.fsu.edu/~jburkardt/examples/mitchell_freefem++/test04b.edp 18 : // 19 : // Licensing: 20 : // 21 : // This code is distributed under the GNU LGPL license. 22 : // 23 : // Modified: 24 : // 25 : // 17 December 2014 26 : // 27 : // Author: 28 : // 29 : // John Burkardt 30 : // 31 : // Reference: 32 : // 33 : // Frederic Hecht, 34 : // Freefem++, 35 : // Third Edition, version 3.22 36 : // 37 : // William Mitchell, 38 : // A collection of 2D elliptic problems for testing adaptive 39 : // grid refinement algorithms, 40 : // Applied Mathematics and Computation, 41 : // Volume 220, 1 September 2013, pages 350-364. 42 : // 43 : border bottom ( t = 0.0, 1.0 ) { x = t; y = 0.0; label = 1; } 44 : border right ( t = 0.0, 1.0 ) { x = 1.0; y = t; label = 1; } 45 : border top ( t = 1.0, 0.0 ) { x = t; y = 1.0; label = 1; } 46 : border left ( t = 1.0, 0.0 ) { x = 0.0; y = t; label = 1; } 47 : // 48 : // Define Th, the triangulation of the "left" side of the boundaries. 49 : // 50 : int n = 20; 51 : mesh Th = buildmesh ( bottom ( n ) + right ( n ) + top ( n ) + left ( n ) ); 52 : // 53 : // Define Vh, the finite element space defined over Th, using P1 basis functions. 54 : // 55 : fespace Vh ( Th, P1 ); 56 : // 57 : // Define U, V, and F, piecewise continuous functions over Th. 58 : // 59 : Vh u; 60 : Vh v; 61 : // 62 : // Define problem parameters. 63 : // 64 : real xc = 0.51; 65 : real yc = 0.117; 66 : real alpha = 100000.0; 67 : // 68 : // Define the right hand side function F. 69 : // 70 : func f = 4.0 * alpha * ( 1.0 - alpha * ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) ) 71 : * exp ( - alpha * ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) ); 72 : // 73 : // Define the boundary function G. 74 : // 75 : func g = exp ( - alpha * ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) ); 76 : // 77 : // Solve the variational problem. 78 : // 79 : solve Laplace ( u, v ) 80 : = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) ) 81 : - int2d ( Th ) ( f * v ) 82 : + on ( 1, u = g ); 83 : // 84 : // Plot the solution. 85 : // 86 : plot ( u, wait = 1, fill = true, ps = "test04b_u.eps" ); 87 : // 88 : // Plot the mesh. 89 : // 90 : plot ( Th, wait = 1, ps = "test04b_mesh.eps" ); 91 : // 92 : // Save the mesh file. 93 : // 94 : savemesh ( Th, "test04b.msh" ); 95 : 96 : sizestack + 1024 =3752 ( 2728 ) -- mesh: Nb of Triangles = 952, Nb of Vertices 517 -- Solve : min -0.012953 max -1.62203e-65 number of required edges : 0 times: compile 0s, execution 0.03s, mpirank:0 CodeAlloc : nb ptr 2697, size :343584 mpirank: 0 Ok: Normal End