load "medit" real L = 4.; //Beam length real W = 1; //Beam height int fixed = 3; //Beam fixed label int free = 2; int loaded = 1; int fya = -2.25e7; int s = 1; //Beam free label border b1(t=-0.025*W,0.025*W){x=0.5*L; y=t; label=loaded;}; border b2(t=0.025*W,0.5*W){x=0.5*L; y=t; label=free;}; border b3(t=0.5*L,-0.5*L){x=t; y=0.5*W; label=free;}; border b4(t=0.5*W,-0.5*W){x=-0.5*L; y=t; label=fixed;}; border b5(t=-0.5*L,0.5*L){x=t; y=-0.5*W; label=free;}; border b6(t=-0.5*W,-0.025*W){x=0.5*L; y=t; label=free;}; mesh Sh = buildmesh(b1(2*s) + b2(24*s) + b3(100*s) + b4(50*s)+ b5(100*s) + b6(24*s)); int a = Sh.nv; //节点数量 //Parameters real Rho = 8000.; //Density real volfrac = 0.5; real chg = 1; real alpha = 0.002; real E1 = 3.e9; //Solid material Young modulus real E0 = E1*1e-9; //empty material Young modulus real Nu = 0.4; //Poisson ratio real Gravity = -9.81; //Gravity int i = 0; //迭代次数 int p = 1; 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)) // 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) ) ) - int1d(Sh,loaded)( fya * vy ) + on(fixed, ux=0, uy=0) ; matrix k = vElasticity(Vh, Vh, solver=sparsesolver); //总刚 real[int] ElasticityBoundary = vElasticity(0, Vh); ux[] = k^-1 * ElasticityBoundary; real[int] u = k^-1 * ElasticityBoundary; real com = [ux,uy]'*k*[ux,uy];