-- 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 : // 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/examples/mitchell_freefem++/test03b.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 : 70 : // 71 : // Set parameters. 72 : // 73 : real nu = 0.3; 74 : real e = 1.0; 75 : real lambda = 0.9085291898461; 76 : real q = -0.2189232362488; 77 : 78 : cout << "\n"; 79 : cout << " Lambda = " << lambda << "\n"; 80 : cout << " Q = " << q << "\n"; 81 : // 82 : // Specify OMEGA, the angle for the slit. 83 : // We can't use omega = 2pi, or even 2pi - 0.05 because the mesher fails. 84 : // 85 : real omega = 8.0 * pi / 4.0 - 0.10; 86 : cout << " Omega = " << omega << "\n"; 87 : // 88 : // Specify P, the relative number of points on each line. 89 : // We can't use P = 3 or the mesher fails. 90 : // 91 : int p = 5; 92 : cout << " Border precision P = " << p << "\n"; 93 : // 94 : // Specify the lines the form the border of the region. 95 : // 96 : border b1 ( t = 0.0, +1.0 ) { x = t; y = 0.0; label = 1; } 97 : real l1 = 1.0; 98 : int n1 = ( round ) ( l1 * p + 1 ); 99 : border b2 ( t = 0.0, +1.0 ) { x = +1.0; y = t; label = 1; } 100 : real l2 = 1.0; 101 : int n2 = ( round ) ( l2 * p + 1 ); 102 : border b3 ( t = 1.0, -1.0 ) { x = t; y = 1.0; label = 1; } 103 : real l3 = 2.0; 104 : int n3 = ( round ) ( l3 * p + 1 ); 105 : border b4 ( t = 1.0, -1.0 ) { x = -1.0; y = t; label = 1; } 106 : real l4 = 2.0; 107 : int n4 = ( round ) ( l4 * p + 1 ); 108 : border b5 ( t = -1.0, +1.0 ) { x = t; y = -1.0; label = 1; } 109 : real l5 = 2.0; 110 : int n5 = ( round ) ( l5 * p + 1 ); 111 : real xc = 1.0; 112 : real yc = tan ( omega ); 113 : real rc = 1.0 / cos ( omega ); 114 : border b6 ( t = -1.0, yc ) { x = 1.0; y = t; label = 1; } 115 : real l6 = yc + 1.0; 116 : int n6 = ( round ) ( l6 * p + 1 ); 117 : border b7 ( t = rc, 0.0 ) { x = t * cos ( omega ); y = t * sin ( omega ); label = 1; } 118 : real l7 = rc; 119 : int n7 = ( round ) ( l7 * p + 1 ); 120 : // 121 : // Create the mesh. 122 : // 123 : mesh Th = buildmesh ( b1 ( n1 ) + b2 ( n2 ) + b3 ( n3 ) + b4 ( n4 ) 124 : + b5 ( n5 ) + b6 ( n6 ) + b7 ( n7 ) ); 125 : // 126 : // Plot the mesh. 127 : // 128 : plot ( Th, wait = 1, ps = "test03b_mesh.eps" ); 129 : // 130 : // Save the mesh file. 131 : // 132 : savemesh ( Th, "test03b.msh" ); 133 : // 134 : // Define Vh, the finite element space defined over Th, using 135 : // a vector of two P1 basis functions. 136 : // 137 : fespace Vh ( Th, [ P1, P1 ] ); 138 : // 139 : // Define U and V, piecewise continuous functions over Th. 140 : // 141 : Vh [ u1, u2 ]; 142 : Vh [ v1, v2 ]; 143 : // 144 : // Define the right hand side function F. 145 : // 146 : func f1 = 0.0; 147 : func f2 = 0.0; 148 : // 149 : // Define the boundary condition function G. 150 : // 151 : // While I would have MUCH PREFERRED to write G1 and G2 as standard C++ 152 : // functions, I could not discover the appropriate way to do so, and so 153 : // I had to pack the whole formula for each into a one-liner. 154 : // 155 : real kappa = 3.0 - 4.0 * nu; 156 : real g = 0.5 * e / ( 1.0 + nu ); 157 : 158 : func g1 = pow ( sqrt ( x * x + y * y ), lambda ) * 159 : ( ( kappa - q * ( lambda + 1.0 ) ) * cos ( lambda * atan2 ( y, x ) ) 160 : - lambda * cos ( ( lambda - 2.0 ) * atan2 ( y, x ) ) ) / 2.0 / g; 161 : 162 : func g2 = pow ( sqrt ( x * x + y * y ), lambda ) * 163 : ( ( kappa - q * ( lambda + 1.0 ) ) * sin ( lambda * atan2 ( y, x ) ) 164 : - lambda * sin ( ( lambda - 2.0 ) * atan2 ( y, x ) ) ) / 2.0 / g; 165 : // 166 : // Set coefficients for the variational problem. 167 : // 168 : real a = + e * ( 1.0 - nu * nu ) / ( 1.0 - 2.0 * nu ); 169 : real b = + e * ( 1.0 - nu * nu ) / ( 2.0 - 2.0 * nu ); 170 : real c = + e * ( 1.0 - nu * nu ) / ( 1.0 - 2.0 * nu ) / ( 2.0 - 2.0 * nu ); 171 : // 172 : // Solve the variational problem. 173 : // u1 = x, u2 = y is accepted as a boundary condition, 174 : // but not u1 = g1, u2 = g2. 175 : // 176 : solve Laplace ( u1, u2, v1, v2 ) 177 : = int2d ( Th ) ( a*dx(u1)*dx(v1) + b*dy(u1)*dy(v1) + c*dx(u2)*dy(v1) ) 178 : + int2d ( Th ) ( b*dx(u2)*dx(v2) + a*dy(u2)*dy(v2) + c*dx(u1)*dy(v2) ) 179 : + on ( 1, u1 = g1, u2 = g2 ); 180 : // 181 : // Plot the solution. 182 : // 183 : plot ( [ u1, u2 ], wait = 1, fill = true, ps = "test03b_u.eps" ); 184 : 185 : sizestack + 1024 =2232 ( 1208 ) 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 times: compile 0s, execution 0.03s, mpirank:0 CodeAlloc : nb ptr 3035, size :354608 mpirank: 0 Ok: Normal End