load "medit" mesh Sh=square(3,3,[8*x,2*y]); int a = Sh.nv; //节点数量 cout << " b = " << " a " << endl; //Parameters real Rho = 8000.; //Density real volfrac = 0.5; real chg = 1; real alpha = 0.002; real E1 = 210.e9; //Solid material Young modulus real E0 = E1*1e-9; //empty material Young modulus real Nu = 0.27; //Poisson ratio real Gravity = -9.81; //Gravity int i = 0; //迭代次数 int p = 3; int Imax = 100; //最大迭代次数 real Mu = 1/(2.*(1.+Nu)); real Lambda = Nu/((1.+ Nu)*(1.-2.*Nu)); fespace Vh(Sh,[P1,P1]); Vh[ux,uy]; fespace Vh1(Sh,P1); Vh1 E,theta,thetaold,thetanew,dttheta,sens,sensn,sensv; int m = 0; for (m = 0; m <= a-1; m++){ theta [][m] = volfrac; } E = theta^p *(E1-E0)+E0; //初始设计 //Macro real sqrt2 = sqrt(2.); macro Epsilon(ux, uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt2] macro Divergence(ux, uy) (dx(ux) + dy(uy)) while (chg > 1e-3 && i <= Imax) { i = i+1;thetaold = theta; varf vElasticity ([ux,uy], [vx, vy]) = int2d(Sh)( E * Lambda * Divergence(ux, uy) * Divergence(vx, vy) + 2. * E * Mu * ( Epsilon(vx, vy)' * Epsilon(ux, uy) ) ) + int2d(Sh)( Rho * Gravity * vy ) + on(4, ux=0, uy=0) ; //求解位移 matrix k = vElasticity(Vh, Vh, solver=sparsesolver); sens = -p*theta^(p-1) * (E1-E0)*(ux,uy)' * k * (ux,uy); solve smoothing(sensn,sensv) = int2d(th)(alpha*( dx(sensn)*dx(sensv)+dy(sensn)*dy(sensv)) + sensn * sensv) - int2d(th) (sens*sensv); //平滑灵敏度 real l1 = 0, real l2 = 10000; real move = 0.1; while ((l2-l1)/(l2+l1) >(1e-4) && l2> 1e-40) { real lmid = 0.5*(l2+l1); thetanew = max(0.,max(theta - move, min(1., min(theta + move, theta*(max(1e-10, -sensn/lmid))^0.3)))); if ((int2d(sh)(thetanew) - int2d(th)(volfrac)) > 0 { l1 = lmid;} else { l2 = lmid;} } theta = thetanew; E = theta ^p *(E1-E0)+E0; real vol = int2d(th)(theta)/int2d(th)(1); chg = dttheta[].max; plot(theta, fill = 1); cout << "iter = "<< i <<";vol = "<