-- FreeFem++ v4.14 (mer. 06 mars 2024 16:59:04 CET - git v4.14-1-g2b2052ae) file : convergence.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // convergence.edp 2 : // 3 : // Discussion: 4 : // 5 : // Solve a simple boundary value problem on the unit square. 6 : // Repeat the process after halving h. 7 : // Report the errors and error convergence rates. 8 : // 9 : // Licensing: 10 : // 11 : // This code is distributed under the MIT license. 12 : // 13 : // Modified: 14 : // 15 : // 11 March 2021 16 : // 17 : // Author: 18 : // 19 : // John Burkardt 20 : // 21 : cout << endl; 22 : cout << "convergence():" << endl; 23 : cout << " FreeFem++ version." << endl; 24 : cout << " Solve a boundary value problem in a square." << endl; 25 : cout << " Repeatedly refine the mesh, and report the L2 e ... : rror norm." << endl; 26 : // 27 : // Declare some convergence variables here. 28 : // 29 : real a; 30 : real aold; 31 : real arate; 32 : real b; 33 : real bold; 34 : real brate; 35 : real c; 36 : real cold; 37 : real crate; 38 : real h; 39 : real hold; 40 : // 41 : // Loop for h = 1/8, 1/16, ... 42 : // 43 : for ( int log2 = 3; log2 <= 7; log2++ ) 44 : { 45 : // 46 : // n: number of elements in each spatial direction. 47 : // 48 : int n = pow ( 2.0, log2 ); 49 : // 50 : // h: characteristic length of each element. 51 : // 52 : if ( 3 < log2 ) 53 : { 54 : hold = h; 55 : } 56 : h = 1.0 / n; 57 : // 58 : // Th: the triangulation of the region. 59 : // 60 : mesh Th = square ( n, n ); 61 : // 62 : // Plot the mesh the first time. 63 : // 64 : if ( log2 == 3 ) 65 : { 66 : plot ( Th, wait = false, 67 : cmm = "Initial mesh of the region", 68 : ps = "convergence_mesh.ps" ); 69 : } 70 : // 71 : // Define Vh, the finite element space. "P1" = piecewise linears. 72 : // 73 : fespace Vh ( Th, P1 ); 74 : // 75 : // ch: the finite element function that will approximate the solution. 76 : // 77 : Vh ch; 78 : // 79 : // vh: the test functions used to project the error. 80 : // 81 : Vh vh; 82 : // 83 : // ce: the exact solution, projected into the finite element space. 84 : // 85 : Vh ce = ( x - x^2 ) * sin ( pi * y ); 86 : // 87 : // f: the right hand side. 88 : // 89 : func f = 90 : ( 1.0 - 2.0 * x ) * sin ( pi * y ) 91 : + ( x - x^2 ) * pi * cos ( pi * y ) 92 : + 2.0 * sin ( pi * y ) 93 : + ( x - x^2 ) * pi^2 * sin ( pi * y ); 94 : // 95 : // Define the problem. 96 : // 97 : problem bvp ( ch, vh ) 98 : = int2d ( Th ) 99 : ( 100 : ( dx(ch) + dy(ch) ) * vh 101 : + dx(ch) * dx(vh) + dy(ch) * dy(vh) 102 : ) 103 : + int2d ( Th ) 104 : ( 105 : - f * vh 106 : ) 107 : + on ( 1, 2, 3, 4, ch = ce ); 108 : // 109 : // Solve the problem. 110 : // 111 : bvp; 112 : // 113 : // Compute error and convergence rates. 114 : // 115 : if ( 3 < log2 ) 116 : { 117 : aold = a; 118 : bold = b; 119 : cold = c; 120 : } 121 : 122 : a = sqrt 123 : ( 124 : int2d ( Th ) 125 : ( 126 : ( ch - ce )^2 127 : ) 128 : ); 129 : 130 : b = sqrt 131 : ( 132 : int2d ( Th ) ( ( dx(ch) - dx(ce) )^2 ) 133 : + int2d ( Th ) ( ( dy(ch) - dy(ce) )^2 ) 134 : ); 135 : 136 : c = sqrt ( a^2 + b^2 ); 137 : 138 : if ( 3 < log2 ) 139 : { 140 : arate = ( log ( aold ) - log ( a ) ) / ( log ( hold ) - log ( h ) ); 141 : brate = ( log ( bold ) - log ( b ) ) / ( log ( hold ) - log ( h ) ); 142 : crate = ( log ( cold ) - log ( c ) ) / ( log ( hold ) - log ( h ) ); 143 : } 144 : 145 : cout << endl; 146 : cout << " H = " << h << endl; 147 : if ( 3 == log2 ) 148 : { 149 : cout << " L2 = " << a << endl; 150 : cout << " H0 = " << b << endl; 151 : cout << " H1 = " << c << endl; 152 : } 153 : else 154 : { 155 : cout << " L2 = " << a << " L2rate = " << arate << endl; 156 : cout << " H0 = " << b << " H0rate = " << brate << endl; 157 : cout << " H1 = " << c << " H1rate = " << crate << endl; 158 : } 159 : } 160 : // 161 : // Terminate. 162 : // 163 : cout << "\n"; 164 : cout << "convergence():\n"; 165 : cout << " Normal end of execution.\n"; 166 : 167 : sizestack + 1024 =2424 ( 1400 ) convergence(): FreeFem++ version. Solve a boundary value problem in a square. Repeatedly refine the mesh, and report the L2 error norm. -- Square mesh : nb vertices =81 , nb triangles = 128 , nb boundary edges 32 -- Solve : min -1.74611e-33 max 0.247077 H = 0.125 L2 = 0.00158174 H0 = 0.00768074 H1 = 0.00784191 -- Square mesh : nb vertices =289 , nb triangles = 512 , nb boundary edges 64 -- Solve : min -2.3862e-34 max 0.249269 H = 0.0625 L2 = 0.000406335 L2rate = 1.96077 H0 = 0.00195573 H0rate = 1.97353 H1 = 0.0019975 H1rate = 1.97301 -- Square mesh : nb vertices =1089 , nb triangles = 2048 , nb boundary edges 128 -- Solve : min -3.09383e-35 max 0.249817 H = 0.03125 L2 = 0.000102284 L2rate = 1.99009 H0 = 0.000491209 H0rate = 1.9933 H1 = 0.000501745 H1rate = 1.99317 -- Square mesh : nb vertices =4225 , nb triangles = 8192 , nb boundary edges 256 -- Solve : min -3.93211e-36 max 0.249954 H = 0.015625 L2 = 2.5615e-05 L2rate = 1.99752 H0 = 0.000122946 H0rate = 1.99832 H1 = 0.000125586 H1rate = 1.99828 -- Square mesh : nb vertices =16641 , nb triangles = 32768 , nb boundary edges 512 -- Solve : min -4.95449e-37 max 0.249989 H = 0.0078125 L2 = 6.4065e-06 L2rate = 1.99938 H0 = 3.07454e-05 H0rate = 1.99958 H1 = 3.14057e-05 H1rate = 1.99957 convergence(): Normal end of execution. times: compile 0.067838s, execution 3.50441s, mpirank:0 CodeAlloc : nb ptr 4165, size :528624 mpirank: 0 Ok: Normal End