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