#include "Sundance.h" /* This code computes sensitivities for the membrane bi-material * optimal design problem that I described in class. The BVP is the Poisson * equation in a square with homogeneous Dirichlet boundary conditions. * The mebrane has a material property of k1 outside a circle, and k2 inside * the circle. The design parameters are the (x,y) coordinates of the * circle, and its radius. The objective function is a weighted sum of the * flexibility of the membrane (measured by the strain energy), and the * "cost" of the material used. The code below computes the sensitivity of * the state variable, i.e. the transverse displacement of the membrane, to * the three design parameters, as well as the sensitivity of the objective * function to those parameters, using finite element approximation of the * appropriate variational forms. A level set method is used to formulate * the problem in a form suitable for a regular mesh. */ int main(int argc, void** argv) { try { TSFCommandLine::verbose() = true; Sundance::init(&argc, &argv); /* create a simple mesh on the rectangle */ int nx = 32; int ny = 32; int npx = 1; int npy = 1; int fill=1; int overlap=1; int femDegree = 2; // degree of finite element approximation int quadOrder = 4; // order of quadrature formula double tol = 1.0e-12; // iterative solver tolerance double alpha = 100.0; // smoothness parameter for heaviside function double beta = 0.01; // regularization parameter double x0 = 0.0; // x coordinate of center of circle of softer material double y0 = 0.0; // y coordinate of center of circle of softer material double r = 1.0; // r radius of circle of softer material double k1 = 5.0; // stiffness of stiffer material double k2 = 1.0; // stiffness of softer material int maxiter = 2000; TSFCommandLine::findDouble("-alpha", alpha); TSFCommandLine::findDouble("-beta", beta); TSFCommandLine::findDouble("-x0", x0); TSFCommandLine::findDouble("-y0", y0); TSFCommandLine::findDouble("-r", r); TSFCommandLine::findDouble("-k1", k1); TSFCommandLine::findDouble("-k2", k2); TSFCommandLine::findInt("-npx", npx); TSFCommandLine::findInt("-npy", npy); TSFCommandLine::findInt("-nx", nx); TSFCommandLine::findInt("-ny", ny); TSFCommandLine::findInt("-maxiter", maxiter); TSFCommandLine::findInt("-fill", fill); TSFCommandLine::findInt("-overlap", overlap); TSFCommandLine::findInt("-femDegree", femDegree); TSFCommandLine::findInt("-quadOrder", quadOrder); TSFOut::printf("npx=%d npy=%d\n", npx, npy); bool useAztec = TSFCommandLine::find("-aztec"); TSFCommandLine::findDouble("-tol", tol); MeshGenerator mesher = new PartitionedRectangleMesher(-2.0, 2.0, npx*nx, npx, -2.0, 2.0, npy*ny, npy); Mesh mesh = mesher.getMesh(); /* define coordinate functions for x and y coordinates */ Expr x = new CoordExpr(0); Expr y = new CoordExpr(1); /* define boundary of rectangle */ CellSet boundary = new BoundaryCellSet(); /* specify the low-level linear algebra representation */ TSFVectorType petra = new PetraVectorType(); /* create a discrete space on the mesh */ TSFVectorSpace discreteSpace = new SundanceVectorSpace(mesh, new Lagrange(femDegree), petra); /* We'll use a discrete function to represent the * source term, providing a test * of our ability to evaluate discrete functions on maximal cells */ Expr f = new DiscreteFunction(discreteSpace, 1.0); /* The material property k(x,y) is defined in terms of a level set function phi, * which gives a circle of softer material (k=1) within a stiffer material (k=5). * The circle has radius r, and is centered at (x0,y0) */ Expr phi = (x-x0)*(x-x0) + (y-y0)*(y-y0) - r*r; Expr chi = 0.5*(1+tanh(alpha*phi)); Expr k = k1*chi + k2*(1-chi); Expr DkDx0 = alpha*(k2-k1)*(x-x0)/pow(cosh(alpha*phi),2); Expr DkDy0 = alpha*(k2-k1)*(y-y0)/pow(cosh(alpha*phi),2); Expr DkDr = alpha*(k2-k1)*r/pow(cosh(alpha*phi),2); /* create symbolic objects for test and unknown functions */ Expr uHat = new TestFunction(new Lagrange(femDegree)); Expr u = new UnknownFunction(new Lagrange(femDegree)); Expr uSensX0 = new UnknownFunction(new Lagrange(femDegree)); Expr uSensY0 = new UnknownFunction(new Lagrange(femDegree)); Expr uSensR = new UnknownFunction(new Lagrange(femDegree)); /* create symbolic differential operators */ Expr dx = new Derivative(0,1); Expr dy = new Derivative(1,1); Expr grad = List(dx, dy); /* Write symbolic weak equation and Neumann and Robin BCs */ Expr weakState = Integral(k*((grad*uHat)*(grad*u)) - f*uHat, new GaussianQuadrature(quadOrder)); /* Write essential BC, u=0 on boundary */ EssentialBC bcState = EssentialBC(boundary, uHat*u); /* Assemble everything into a problem object, with a specification that * Petra be used as the low-level linear algebra representation */ StaticLinearProblem stateProb(mesh, weakState, bcState, uHat, u, petra); /* create a preconditioner and solver */ TSFHashtable azOptions; TSFHashtable azParams; TSFLinearSolver solver; if (useAztec) { azOptions.put(AZ_solver, AZ_gmres); azOptions.put(AZ_precond, AZ_dom_decomp); azOptions.put(AZ_subdomain_solve, AZ_ilu); azOptions.put(AZ_graph_fill, fill); azParams.put(AZ_tol, tol); azOptions.put(AZ_max_iter, maxiter); solver = new AZTECSolver(azOptions, azParams); } else { TSFPreconditionerFactory precond = new ILUKPreconditionerFactory(fill, overlap); solver = new BICGSTABSolver(precond, tol, maxiter); } /* solve the problem and return the solution as a symbolic object */ Expr state = stateProb.solve(solver); // Expr exactSoln = 0.5*x*x + y/3.0; /* write to matlab */ FieldWriter writer = new MatlabWriter("state.dat"); writer.writeField(state); /* solve sensitivity equation for X0 */ Expr weakSensX0 = Integral(k*((grad*uHat)*(grad*uSensX0)) + DkDx0*((grad*state)*(grad*uHat)), new GaussianQuadrature(quadOrder)); EssentialBC bcSensX0 = EssentialBC(boundary, uHat*uSensX0); StaticLinearProblem sensX0Prob(mesh, weakSensX0, bcSensX0, uHat, uSensX0, petra); Expr sensX0 = sensX0Prob.solve(solver); FieldWriter writer2 = new MatlabWriter("sensX0.dat"); writer2.writeField(sensX0); /* solve sensitivity equation for Y0 */ Expr weakSensY0 = Integral(k*((grad*uHat)*(grad*uSensY0)) + DkDy0*((grad*state)*(grad*uHat)), new GaussianQuadrature(quadOrder)); EssentialBC bcSensY0 = EssentialBC(boundary, uHat*uSensY0); StaticLinearProblem sensY0Prob(mesh, weakSensY0, bcSensY0, uHat, uSensY0, petra); Expr sensY0 = sensY0Prob.solve(solver); FieldWriter writer3 = new MatlabWriter("sensY0.dat"); writer3.writeField(sensY0); /* solve sensitivity equation for R */ Expr weakSensR = Integral(k*((grad*uHat)*(grad*uSensR)) + DkDr*((grad*state)*(grad*uHat)), new GaussianQuadrature(quadOrder)); EssentialBC bcSensR = EssentialBC(boundary, uHat*uSensR); StaticLinearProblem sensRProb(mesh, weakSensR, bcSensR, uHat, uSensR, petra); Expr sensR = sensRProb.solve(solver); FieldWriter writer4 = new MatlabWriter("sensR.dat"); writer4.writeField(sensR); // Compute objective function Expr strainEnergy = 0.5*(grad*state)*(grad*state); Expr cost = beta*k; double objective = strainEnergy.integral(mesh) + cost.integral(mesh); cout << " strainEnergy = " << strainEnergy.integral(mesh) << " cost = " << cost.integral(mesh) << " objective = " << objective << endl; // Compute objective function gradient w.r.t. X0 Expr strainEnergyGradX0 = (grad*state)*(grad*sensX0); Expr costGradX0 = beta*DkDx0; double objectiveGradX0 = strainEnergyGradX0.integral(mesh) + costGradX0.integral(mesh); cout << " strainEnergyGradX0 = " << strainEnergyGradX0.integral(mesh) << " costGradX0 = " << costGradX0.integral(mesh) << " objectiveGradX0 = " << objectiveGradX0 << endl; // Compute objective function gradient w.r.t. Y0 Expr strainEnergyGradY0 = (grad*state)*(grad*sensY0); Expr costGradY0 = beta*DkDy0; double objectiveGradY0 = strainEnergyGradY0.integral(mesh) + costGradY0.integral(mesh); cout << " strainEnergyGradY0 = " << strainEnergyGradY0.integral(mesh) << " costGradY0 = " << costGradY0.integral(mesh) << " objectiveGradY0 = " << objectiveGradY0 << endl; // Compute objective function gradient w.r.t. R Expr strainEnergyGradR = (grad*state)*(grad*sensR); Expr costGradR = beta*DkDr; double objectiveGradR = strainEnergyGradR.integral(mesh) + costGradR.integral(mesh); cout << " strainEnergyGradR = " << strainEnergyGradR.integral(mesh) << " costGradR = " << costGradR.integral(mesh) << " objectiveGradR = " << objectiveGradR << endl; } catch(exception& e) { Sundance::handleError(e, __FILE__); } Sundance::finalize(); }