#include "Sundance.h" /* * This is code for solving the minimal surface problem, based * on code written by Clemens Kadow. The problem is to minimize the * surface area = integral((grad u . grad u + 1)^0.5) over Omega, * subject to u = g on Gamma. */ int main(int argc, void** argv) { try { Sundance::init(argc, argv); int nx = 32; // number of elements in x direction int ny = 32; // number of elements in y direction int npx = 1; // number of processors in x direction int npy = 1; // number of processors in y direction int fill=1; // parameter that goes into preconditioner // larger fill means fewer iterations, but more work // per iteration; sweet spot tends to be ~ 1 or 2 int overlap=1; // only meaningful for parallel solver double ax = -1.0; // left x coordinate defining domain double ay = -1.0; // right x coordinate defining domain double bx = 1.0; // bottom y coordinate defining domain double by = 1.0; // top y coordinate defining domain double delta =1.0e-2; int period = 10; // this is the period of a sinusoidal function that // is the initial guess for the shape of the surface, // and is also the boundary condition int order = 1; double lintol = 1.0e-10; // linear solver tolerance double ftol = 1.0e-10; // newton solver tolerance on residual double stol = 1.0e-10; // newton solver tolerance on change in resi int maxiterlin = 100; // max # of linear iterations int maxiter = 100; // max # of nonlinear iterations // the following lines let you input any of the listed parameters // via command line swtiches TSFCommandLine::findInt("-nx", nx); TSFCommandLine::findInt("-ny", ny); TSFCommandLine::findInt("-npx", npx); TSFCommandLine::findInt("-npy", npy); TSFCommandLine::findInt("-maxiter", maxiter); TSFCommandLine::findInt("-maxiterlin",maxiterlin); TSFCommandLine::findInt("-fill", fill); TSFCommandLine::findInt("-overlap", overlap); TSFCommandLine::findDouble("-lintol",lintol); TSFCommandLine::findDouble("-ftol",ftol); TSFCommandLine::findDouble("-stol",stol); TSFCommandLine::findDouble("-delta",delta); TSFCommandLine::findInt("-period",period); TSFCommandLine::findInt("-order",order); /* Create a mesh object. In this example, we will use a built-in method to create a uniform mesh on the unit line. For non-rectangular geomtery, we would use a mesher to create a mesh, and then read the mesh using a MeshReader object. */ MeshGenerator mesher = new PartitionedRectangleMesher(ax, bx, nx, npx, ay, by, ny, npy); //MeshReader mesher = new ShewchukMeshReader("BalticSea/balticsea30"); Mesh mesh = mesher.getMesh();//.getSubmesh(); //cout << "Read mesh from file." << endl; /* Define a symbolic object to represent the coordinate functions. */ Expr x = new CoordExpr(0); Expr y = new CoordExpr(1); /* * Define a cell set that contains all boundary cells. * For the minimal surface problem, we'll just use a single * expression for our BC, so we don't need to distinguish * the four different boundaries (top, bottom, left, right). * Uncomment lines below if you want different BCs on the the four * boundaries. */ CellSet boundary = new BoundaryCellSet(); cout << "Specified boundary" << endl; /* Define a cell sets of all four boundaries */ //CellSet left = boundary.subset( fabs(x - ax) < 1.0e-10 ); //CellSet right = boundary.subset( fabs(x - bx) < 1.0e-10 ); //CellSet top = boundary.subset( fabs(y - by) < 1.0e-10 ); //CellSet bottom = boundary.subset( fabs(y - ay) < 1.0e-10 ); /* Generate a discrete vector space of first-order * Lagrange polynomials. You can change to second-order if you wish. */ TSFVectorSpace discreteSpace = new SundanceVectorSpace(mesh,new Lagrange(order)); cout << "Generated discrete Spaces" << endl; /* Define an unknown function and its variation. The constructor argument is the basis family with which the function will be represented, in this case first-order Lagrange (nodal) polynomials. */ //Expr initial = 1/(x+delta) * 1/(y+delta); //Expr initial = -(0.1-x)*(0.1+x)*(0.05-x)*(0.05+x)*(0.1-y) //+ (0.1-x)*(0.1+x)*(x)*(0.1+y) //+ (0.1-y)*(y)*(0.1+y)*(0.1-x) //+ (0.1-y)*(y)*(0.1+y)*(0.1+x); Expr initial = cos(period*x)*cos(period*y); /* basis, initial guess, and test function for the displacement */ Expr u = new UnknownFunction(new Lagrange(order)); cout << "Generated u" << endl; // Expr u0 = DiscreteFunction::discretize(discreteSpace, initial ); Expr u0 = new DiscreteFunction(discreteSpace, initial); cout << "Generated u0" << endl; Expr v = u.variation(); Expr bcs = u0; cout << "Generated bcs" << endl; /* Define differentiation operators */ Expr dx = new Derivative(0); Expr dy = new Derivative(1); /* the gradient operator applied to a scalar */ Expr scalGrad = List(dx,dy); Expr one = new ParameterExpr(1.0); /* Objective functional -- the surface area */ Expr F = Integral( sqrt(one+(scalGrad*u)*(scalGrad*u)) ); //cout << "surface area functional: " << F << endl; // Now let's take the variation of the objective w.r.t. u // This produces the "nonLinEqn" expression // Expr nonLinEqn = F.variation(u); // or we could have just done this by hand Expr nonLinEqn = Integral(1/sqrt((scalGrad*u) * (scalGrad*u) + one ) * (scalGrad*u)*(scalGrad*v) ); //cout << "eqn = " << nonLinEqn << endl; /* We want to solve nonLinEqn = 0 * So we tell sundance to linearize the equation w.r.t u, * about the function u0, which produces linEqn */ Expr linEqn = nonLinEqn.linearization(u,u0); // Again, we could have done this linearization by hand // (and would have to, if we want an approximate Hessian) //cout << linEqn << endl; /* In doing the linearization, sundance produced the * differential of u, which we will call du. Note that this is * Newton search direction */ Expr du = u.differential(); /* * Now specify the weak form of the boundary conditions on the * "boundary" (this would be top, bottom, left, right, if you had * distinguished above. */ EssentialBC nonLinBC = EssentialBC(boundary, (u-bcs)*v); // && EssentialBC(right, (u-bcs)*v) //&& EssentialBC(top, (u-bcs)*v) //&& EssentialBC(bottom, (u-bcs)*v); /* Even though we have linear BCs, they are in terms of u, and * they still need to be linearized, i.e. put into a form that * involves the differential du. This is accomplished by the * following line, where we are linearizing the variable u about * function u0. This produces the linBC expression. */ EssentialBC linBC = nonLinBC.linearization(u,u0); // Or else we could directly create te linearized BCs ourselves; just be careful to distinguish // the current iterate u0 and the Newton step du, as follows: //EssentialBC linBC = EssentialBC(boundary, du*v + (u0-bcs)*v); /* Create a solver object: stablized biconjugate gradient solver */ TSFPreconditionerFactory prec = new ILUKPreconditionerFactory(fill); TSFLinearSolver solver = new BICGSTABSolver(prec, lintol, maxiterlin); /* * Combine the geometry, the variational form, the BCs, and the solver * to form a complete problem. First we create the linear problem * (with mesh, linearized equation, linearized boundary conditions, * test function, unknoiwn function), and then the nonlinear problem * (which includes the just-defined linear problem, in addition to * initial guess and linear solver). */ StaticLinearProblem linProb(mesh, linEqn, linBC, v, du); NewtonLinearization nonLinProb(linProb, u0, solver); // Next we write out the initial guess in a form that the matlab code // plot2D can read char fInitName[100]; sprintf(fInitName, "ms-hInit%d.dat.%d", period, MPIComm::world().getRank()); FieldWriter initWriter = new MatlabWriter(fInitName); initWriter.writeField(u0); /* now we actually solve the nonlinear problem, using the nonlinear * the nonlinear problem (nonlinprob) */ Expr soln = nonLinProb.solve(NewtonSolver(solver, maxiter, ftol, stol)); // finally write out the solution in matlab form char fName[100]; sprintf(fName, "ms-h%d.dat.%d", period, MPIComm::world().getRank()); FieldWriter writer = new MatlabWriter(fName); writer.writeField(u0); } catch(exception& e) { TSFOut::println(e.what()); Testing::crash(__FILE__); Testing::timeStamp(__FILE__, __DATE__, __TIME__); } Sundance::finalize(); }