-- FreeFem++ v4.14 (mer. 06 mars 2024 16:59:04 CET - git v4.14-1-g2b2052ae) file : mitchell_03b.edp Load: lg_fem lg_mesh lg_mesh3 eigenvalue 1 : // mitchell_03b.edp 2 : // 3 : // Discussion: 4 : // 5 : // Coupled elasticity equations on a square with a slit. 6 : // 7 : // a * u1xx + b * u1yy + c * u2xy = f1 8 : // b * u2xx + a * u2yy + c * u1xy = f2 9 : // 10 : // f1 = 0 11 : // f2 = 0 12 : // u1 = g1 on boundary 13 : // u2 = g2 on boundary 14 : // 15 : // a = - E * ( 1 - nu^2 ) / ( 1 - 2 * nu ) 16 : // b = - E * ( 1 - nu^2 ) / ( 2 - 2 * nu ) 17 : // c = - E * ( 1 - nu^2 ) / ( 1 - 2 * nu ) / ( 2 - 2 * nu ) 18 : // 19 : // E = 1 20 : // nu = 0.3 21 : // 22 : // kappa = 3 - 4 * nu 23 : // g = 0.5 * e / ( 1 + nu ) 24 : // 25 : // c1 = ( kappa - q * ( lambda + 1 ) * cos ( lambda * theta ) 26 : // c2 = lambda * cos ( ( lambda - 2.0 ) * theta ) 27 : // g1(r,theta) = r^lambda * ( c1 - c2 ) / 2 / g 28 : // s1 = ( kappa - q * ( lambda + 1 ) * sin ( lambda * theta ) 29 : // s2 = lambda * sin ( ( lambda - 2.0 ) * theta ) 30 : // g2(r,theta) = r^lambda * ( s1 - s2 ) / 2 / g 31 : // 32 : // Suggested Parameter Values: 33 : // 34 : // lambda = 0.5444837367825, q = 0.5430755788367 35 : // lambda = 0.9085291898461, q = -0.2189232362488 36 : // 37 : // Licensing: 38 : // 39 : // This code is distributed under the MIT license. 40 : // 41 : // Modified: 42 : // 43 : // 15 February 2015 44 : // 45 : // Author: 46 : // 47 : // John Burkardt 48 : // 49 : // Reference: 50 : // 51 : // Mark Ainsworth, Bill Senior, 52 : // Aspects of an adaptive hp-finite element method: 53 : // Adaptive strategy for conforming approximation and efficient solvers, 54 : // Computer Methods in Applied Mechanics and Engineering, 55 : // Volume 150, 1997, pages 65-87, 56 : // 57 : // Frederic Hecht, 58 : // Freefem++, 59 : // Third Edition, version 3.22 60 : // 61 : // William Mitchell, 62 : // A collection of 2D elliptic problems for testing adaptive 63 : // grid refinement algorithms, 64 : // Applied Mathematics and Computation, 65 : // Volume 220, 1 September 2013, pages 350-364. 66 : // 67 : cout << "\n"; 68 : cout << "mitchell_03b():\n"; 69 : cout << " FreeFem++ version\n"; 70 : cout << " Mitchell linear elasticity problem.\n"; 71 : // 72 : // Set parameters. 73 : // 74 : real nu = 0.3; 75 : real e = 1.0; 76 : real lambda = 0.9085291898461; 77 : real q = -0.2189232362488; 78 : 79 : cout << "\n"; 80 : cout << " Lambda = " << lambda << "\n"; 81 : cout << " Q = " << q << "\n"; 82 : // 83 : // Specify OMEGA, the angle for the slit. 84 : // We can't use omega = 2pi, or even 2pi - 0.05 because the mesher fails. 85 : // 86 : real omega = 8.0 * pi / 4.0 - 0.10; 87 : cout << " Omega = " << omega << "\n"; 88 : // 89 : // Specify P, the relative number of points on each line. 90 : // We can't use P = 3 or the mesher fails. 91 : // 92 : int p = 5; 93 : cout << " Border precision P = " << p << "\n"; 94 : // 95 : // Specify the lines the form the border of the region. 96 : // 97 : border b1 ( t = 0.0, +1.0 ) { x = t; y = 0.0; label = 1; } 98 : real l1 = 1.0; 99 : int n1 = ( round ) ( l1 * p + 1 ); 100 : border b2 ( t = 0.0, +1.0 ) { x = +1.0; y = t; label = 1; } 101 : real l2 = 1.0; 102 : int n2 = ( round ) ( l2 * p + 1 ); 103 : border b3 ( t = 1.0, -1.0 ) { x = t; y = 1.0; label = 1; } 104 : real l3 = 2.0; 105 : int n3 = ( round ) ( l3 * p + 1 ); 106 : border b4 ( t = 1.0, -1.0 ) { x = -1.0; y = t; label = 1; } 107 : real l4 = 2.0; 108 : int n4 = ( round ) ( l4 * p + 1 ); 109 : border b5 ( t = -1.0, +1.0 ) { x = t; y = -1.0; label = 1; } 110 : real l5 = 2.0; 111 : int n5 = ( round ) ( l5 * p + 1 ); 112 : real xc = 1.0; 113 : real yc = tan ( omega ); 114 : real rc = 1.0 / cos ( omega ); 115 : border b6 ( t = -1.0, yc ) { x = 1.0; y = t; label = 1; } 116 : real l6 = yc + 1.0; 117 : int n6 = ( round ) ( l6 * p + 1 ); 118 : border b7 ( t = rc, 0.0 ) { x = t * cos ( omega ); y = t * sin ( omega ); label = 1; } 119 : real l7 = rc; 120 : int n7 = ( round ) ( l7 * p + 1 ); 121 : // 122 : // Create the mesh. 123 : // 124 : mesh Th = buildmesh ( b1 ( n1 ) + b2 ( n2 ) + b3 ( n3 ) + b4 ( n4 ) 125 : + b5 ( n5 ) + b6 ( n6 ) + b7 ( n7 ) ); 126 : // 127 : // Plot the mesh. 128 : // 129 : plot ( Th, wait = true, ps = "mitchell_03b_mesh.eps" ); 130 : // 131 : // Save the mesh file. 132 : // 133 : savemesh ( Th, "mitchell_03b.msh" ); 134 : // 135 : // Define Vh, the finite element space defined over Th, using 136 : // a vector of two P1 basis functions. 137 : // 138 : fespace Vh ( Th, [ P1, P1 ] ); 139 : // 140 : // Define U and V, piecewise continuous functions over Th. 141 : // 142 : Vh [ u1, u2 ]; 143 : Vh [ v1, v2 ]; 144 : // 145 : // Define the right hand side function F. 146 : // 147 : func f1 = 0.0; 148 : func f2 = 0.0; 149 : // 150 : // Define the boundary condition function G. 151 : // 152 : // While I would have MUCH PREFERRED to write G1 and G2 as standard C++ 153 : // functions, I could not discover the appropriate way to do so, and so 154 : // I had to pack the whole formula for each into a one-liner. 155 : // 156 : real kappa = 3.0 - 4.0 * nu; 157 : real g = 0.5 * e / ( 1.0 + nu ); 158 : 159 : func g1 = pow ( sqrt ( x * x + y * y ), lambda ) * 160 : ( ( kappa - q * ( lambda + 1.0 ) ) * cos ( lambda * atan2 ( y, x ) ) 161 : - lambda * cos ( ( lambda - 2.0 ) * atan2 ( y, x ) ) ) / 2.0 / g; 162 : 163 : func g2 = pow ( sqrt ( x * x + y * y ), lambda ) * 164 : ( ( kappa - q * ( lambda + 1.0 ) ) * sin ( lambda * atan2 ( y, x ) ) 165 : - lambda * sin ( ( lambda - 2.0 ) * atan2 ( y, x ) ) ) / 2.0 / g; 166 : // 167 : // Set coefficients for the variational problem. 168 : // 169 : real a = + e * ( 1.0 - nu * nu ) / ( 1.0 - 2.0 * nu ); 170 : real b = + e * ( 1.0 - nu * nu ) / ( 2.0 - 2.0 * nu ); 171 : real c = + e * ( 1.0 - nu * nu ) / ( 1.0 - 2.0 * nu ) / ( 2.0 - 2.0 * nu ); 172 : // 173 : // Solve the variational problem. 174 : // u1 = x, u2 = y is accepted as a boundary condition, 175 : // but not u1 = g1, u2 = g2. 176 : // 177 : solve Laplace ( u1, u2, v1, v2 ) 178 : = int2d ( Th ) ( a*dx(u1)*dx(v1) + b*dy(u1)*dy(v1) + c*dx(u2)*dy(v1) ) 179 : + int2d ( Th ) ( b*dx(u2)*dx(v2) + a*dy(u2)*dy(v2) + c*dx(u1)*dy(v2) ) 180 : + on ( 1, u1 = g1, u2 = g2 ); 181 : // 182 : // Plot the solution. 183 : // 184 : plot ( [ u1, u2 ], wait = true, fill = true, ps = "mitchell_03b_u.eps" ); 185 : // 186 : // Terminate. 187 : // 188 : cout << "\n"; 189 : cout << "mitchell_03b():\n"; 190 : cout << " Normal end of execution.\n"; 191 : 192 : exit ( 0 ); 193 : sizestack + 1024 =2232 ( 1208 ) mitchell_03b(): FreeFem++ version Mitchell linear elasticity problem. Lambda = 0.908529 Q = -0.218923 Omega = 6.18319 Border precision P = 5 -- mesh: Nb of Triangles = 292, Nb of Vertices 175 number of required edges : 0 -- Solve : min -4.19899 max 4.19899 mitchell_03b(): Normal end of execution. current line = 192 exit(0) err code 0 , mpirank 0 CodeAlloc : nb ptr 4382, size :538488 mpirank: 0 Ok: Normal End