1 // MFEM Example 9
2 // Nonlinear Constrained Optimization Modification
3 //
4 // Compile with: make ex9
5 //
6 // Sample runs:
7 // ex9 -m ../../data/periodic-segment.mesh -r 3 -p 0 -o 2 -dt 0.002 -opt 1
8 // ex9 -m ../../data/periodic-segment.mesh -r 3 -p 0 -o 2 -dt 0.002 -opt 2
9 //
10 // ex9 -m ../../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10 -opt 1
11 // ex9 -m ../../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10 -opt 2
12 //
13 // ex9 -m ../../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9 -opt 1
14 // ex9 -m ../../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9 -opt 2
15 //
16 // ex9 -m ../../data/amr-quad.mesh -p 1 -r 1 -dt 0.002 -tf 9 -opt 1
17 // ex9 -m ../../data/amr-quad.mesh -p 1 -r 1 -dt 0.002 -tf 9 -opt 2
18 //
19 // ex9 -m ../../data/disc-nurbs.mesh -p 1 -r 2 -dt 0.005 -tf 9 -opt 1
20 // ex9 -m ../../data/disc-nurbs.mesh -p 1 -r 2 -dt 0.005 -tf 9 -opt 2
21 //
22 // ex9 -m ../../data/disc-nurbs.mesh -p 2 -r 2 -dt 0.01 -tf 9 -opt 1
23 // ex9 -m ../../data/disc-nurbs.mesh -p 2 -r 2 -dt 0.01 -tf 9 -opt 2
24 //
25 // ex9 -m ../../data/periodic-square.mesh -p 3 -r 3 -dt 0.0025 -tf 9 -opt 1
26 // ex9 -m ../../data/periodic-square.mesh -p 3 -r 3 -dt 0.0025 -tf 9 -opt 2
27 //
28 // ex9 -m ../../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8 -opt 1
29 // ex9 -m ../../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8 -opt 2
30
31 // Description: This example modifies the standard MFEM ex9 by adding nonlinear
32 // constrained optimization capabilities through the SLBQP and
33 // HIOP solvers. It demonstrates how a user can define a custom
34 // class OptimizationProblem that includes linear/nonlinear
35 // equality/inequality constraints. This optimization is applied
36 // as post-processing to the solution of the transport equation.
37 //
38 // Description of ex9:
39 // This example code solves the time-dependent advection equation
40 // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
41 // u0(x)=u(0,x) is a given initial condition.
42 //
43 // The example demonstrates the use of Discontinuous Galerkin (DG)
44 // bilinear forms in MFEM (face integrators), the use of explicit
45 // ODE time integrators, the definition of periodic boundary
46 // conditions through periodic meshes, as well as the use of GLVis
47 // for persistent visualization of a time-evolving solution. The
48 // saving of time-dependent data files for external visualization
49 // with VisIt (visit.llnl.gov) is also illustrated.
50
51 #include "mfem.hpp"
52 #include <fstream>
53 #include <iostream>
54
55 #ifndef MFEM_USE_HIOP
56 #error This example requires that MFEM is built with MFEM_USE_HIOP=YES
57 #endif
58
59 using namespace std;
60 using namespace mfem;
61
62 // Choice for the problem setup. The fluid velocity, initial condition and
63 // inflow boundary condition are chosen based on this parameter.
64 int problem;
65
66 // Nonlinear optimizer.
67 int optimizer_type;
68
69 // Velocity coefficient
70 bool invert_velocity = false;
71 void velocity_function(const Vector &x, Vector &v);
72
73 // Initial condition
74 double u0_function(const Vector &x);
75
76 // Inflow boundary condition
77 double inflow_function(const Vector &x);
78
79 // Mesh bounding box
80 Vector bb_min, bb_max;
81
82 /// Computes C(x) = sum w_i x_i, where w is a given Vector.
83 class LinearScaleOperator : public Operator
84 {
85 private:
86 const Vector &w;
87 mutable DenseMatrix grad;
88
89 public:
LinearScaleOperator(const Vector & weight)90 LinearScaleOperator(const Vector &weight)
91 : Operator(1, weight.Size()), w(weight), grad(1, width)
92 {
93 for (int i = 0; i < width; i++) { grad(0, i) = w(i); }
94 }
95
Mult(const Vector & x,Vector & y) const96 virtual void Mult(const Vector &x, Vector &y) const
97 {
98 y(0) = w * x;
99 }
100
GetGradient(const Vector & x) const101 virtual Operator &GetGradient(const Vector &x) const
102 {
103 return grad;
104 }
105 };
106
107 /// Nonlinear monotone bounded operator to test nonlinear inequality constraints
108 /// Computes D(x) = tanh(sum(x_i)).
109 class TanhSumOperator : public Operator
110 {
111 private:
112 mutable DenseMatrix grad;
113
114 public:
TanhSumOperator(int size)115 TanhSumOperator(int size) : Operator(1, size), grad(1, width) { }
116
Mult(const Vector & x,Vector & y) const117 virtual void Mult(const Vector &x, Vector &y) const
118 {
119 y(0) = std::tanh(x.Sum());
120 }
121
GetGradient(const Vector & x) const122 virtual Operator &GetGradient(const Vector &x) const
123 {
124 const double ts = std::tanh(x.Sum());
125 const double dtanh = 1.0 - ts * ts;
126 for (int i = 0; i < width; i++) { grad(0, i) = dtanh; }
127 return grad;
128 }
129 };
130
131 /** Monotone and conservative a posteriori correction for transport solutions:
132 * Find x that minimizes 0.5 || x - x_HO ||^2, subject to
133 * sum w_i x_i = mass,
134 * tanh(sum(x_i_min)) <= tanh(sum(x_i)) <= tanh(sum(x_i_max)),
135 * x_i_min <= x_i <= x_i_max,
136 */
137 class OptimizedTransportProblem : public OptimizationProblem
138 {
139 private:
140 const Vector &x_HO;
141 Vector massvec, d_lo, d_hi;
142 const LinearScaleOperator LSoper;
143 const TanhSumOperator TSoper;
144
145 public:
OptimizedTransportProblem(const Vector & xho,const Vector & w,double mass,const Vector & xmin,const Vector & xmax)146 OptimizedTransportProblem(const Vector &xho, const Vector &w, double mass,
147 const Vector &xmin, const Vector &xmax)
148 : OptimizationProblem(xho.Size(), NULL, NULL),
149 x_HO(xho), massvec(1), d_lo(1), d_hi(1),
150 LSoper(w), TSoper(w.Size())
151 {
152 C = &LSoper;
153 massvec(0) = mass;
154 SetEqualityConstraint(massvec);
155
156 D = &TSoper;
157 d_lo(0) = std::tanh(xmin.Sum());
158 d_hi(0) = std::tanh(xmax.Sum());
159 MFEM_ASSERT(d_lo(0) < d_hi(0),
160 "The bounds produce an infeasible optimization problem");
161 SetInequalityConstraint(d_lo, d_hi);
162
163 SetSolutionBounds(xmin, xmax);
164 }
165
CalcObjective(const Vector & x) const166 virtual double CalcObjective(const Vector &x) const
167 {
168 double res = 0.0;
169 for (int i = 0; i < input_size; i++)
170 {
171 const double d = x(i) - x_HO(i);
172 res += d * d;
173 }
174 return 0.5 * res;
175 }
176
CalcObjectiveGrad(const Vector & x,Vector & grad) const177 virtual void CalcObjectiveGrad(const Vector &x, Vector &grad) const
178 {
179 for (int i = 0; i < input_size; i++) { grad(i) = x(i) - x_HO(i); }
180 }
181 };
182
183
184 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
185 form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
186 and advection matrices, and b describes the flow on the boundary. This can
187 be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
188 used to evaluate the right-hand side. */
189 class FE_Evolution : public TimeDependentOperator
190 {
191 private:
192 SparseMatrix &M, &K;
193 const Vector &b;
194 DSmoother M_prec;
195 CGSolver M_solver;
196
197 mutable Vector z;
198
199 double dt;
200 BilinearForm &bf;
201 Vector &M_rowsums;
202
203 public:
204 FE_Evolution(SparseMatrix &M_, SparseMatrix &K_, const Vector &b_,
205 BilinearForm &bf_, Vector &M_rs);
206
SetTimeStep(double dt_)207 void SetTimeStep(double dt_) { dt = dt_; }
SetK(SparseMatrix & K_)208 void SetK(SparseMatrix &K_) { K = K_; }
209 virtual void Mult(const Vector &x, Vector &y) const;
210
~FE_Evolution()211 virtual ~FE_Evolution() { }
212 };
213
214
main(int argc,char * argv[])215 int main(int argc, char *argv[])
216 {
217 // 1. Parse command-line options.
218 problem = 0;
219 optimizer_type = 2;
220 const char *mesh_file = "../../data/periodic-hexagon.mesh";
221 int ref_levels = 2;
222 int order = 3;
223 int ode_solver_type = 3;
224 double t_final = 1.0;
225 double dt = 0.01;
226 bool visualization = true;
227 bool visit = false;
228 bool binary = false;
229 int vis_steps = 5;
230
231 int precision = 8;
232 cout.precision(precision);
233
234 OptionsParser args(argc, argv);
235 args.AddOption(&mesh_file, "-m", "--mesh",
236 "Mesh file to use.");
237 args.AddOption(&problem, "-p", "--problem",
238 "Problem setup to use. See options in velocity_function().");
239 args.AddOption(&ref_levels, "-r", "--refine",
240 "Number of times to refine the mesh uniformly.");
241 args.AddOption(&order, "-o", "--order",
242 "Order (degree) of the finite elements.");
243 args.AddOption(&optimizer_type, "-opt", "--optimizer",
244 "Nonlinear optimizer: 1 - SLBQP,\n\t"
245 " 2 - HIOP.");
246 args.AddOption(&ode_solver_type, "-s", "--ode-solver",
247 "ODE solver: 1 - Forward Euler,\n\t"
248 " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
249 args.AddOption(&t_final, "-tf", "--t-final",
250 "Final time; start time is 0.");
251 args.AddOption(&dt, "-dt", "--time-step",
252 "Time step.");
253 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
254 "--no-visualization",
255 "Enable or disable GLVis visualization.");
256 args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
257 "--no-visit-datafiles",
258 "Save data files for VisIt (visit.llnl.gov) visualization.");
259 args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
260 "--ascii-datafiles",
261 "Use binary (Sidre) or ascii format for VisIt data files.");
262 args.AddOption(&vis_steps, "-vs", "--visualization-steps",
263 "Visualize every n-th timestep.");
264 args.Parse();
265 if (!args.Good())
266 {
267 args.PrintUsage(cout);
268 return 1;
269 }
270 args.PrintOptions(cout);
271
272 // 2. Read the mesh from the given mesh file. We can handle geometrically
273 // periodic meshes in this code.
274 Mesh *mesh = new Mesh(mesh_file, 1, 1);
275 int dim = mesh->Dimension();
276
277 // 3. Define the ODE solver used for time integration. Several explicit
278 // Runge-Kutta methods are available.
279 ODESolver *ode_solver = NULL;
280 switch (ode_solver_type)
281 {
282 case 1: ode_solver = new ForwardEulerSolver; break;
283 case 2: ode_solver = new RK2Solver(1.0); break;
284 case 3: ode_solver = new RK3SSPSolver; break;
285 case 4: ode_solver = new RK4Solver; break;
286 case 6: ode_solver = new RK6Solver; break;
287 default:
288 cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
289 delete mesh;
290 return 3;
291 }
292
293 // 4. Refine the mesh to increase the resolution. In this example we do
294 // 'ref_levels' of uniform refinement, where 'ref_levels' is a
295 // command-line parameter. If the mesh is of NURBS type, we convert it to
296 // a (piecewise-polynomial) high-order mesh.
297 for (int lev = 0; lev < ref_levels; lev++)
298 {
299 mesh->UniformRefinement();
300 }
301 if (mesh->NURBSext)
302 {
303 mesh->SetCurvature(max(order, 1));
304 }
305 mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
306
307 // 5. Define the discontinuous DG finite element space of the given
308 // polynomial order on the refined mesh.
309 DG_FECollection fec(order, dim, BasisType::Positive);
310 FiniteElementSpace fes(mesh, &fec);
311
312 cout << "Number of unknowns: " << fes.GetVSize() << endl;
313
314 // 6. Set up and assemble the bilinear and linear forms corresponding to the
315 // DG discretization. The DGTraceIntegrator involves integrals over mesh
316 // interior faces.
317 VectorFunctionCoefficient velocity(dim, velocity_function);
318 FunctionCoefficient inflow(inflow_function);
319 FunctionCoefficient u0(u0_function);
320
321 BilinearForm m(&fes);
322 m.AddDomainIntegrator(new MassIntegrator);
323 BilinearForm k(&fes);
324 k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
325 k.AddInteriorFaceIntegrator(
326 new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
327 k.AddBdrFaceIntegrator(
328 new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
329
330 LinearForm b(&fes);
331 b.AddBdrFaceIntegrator(
332 new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
333
334 m.Assemble();
335 m.Finalize();
336 int skip_zeros = 0;
337 k.Assemble(skip_zeros);
338 k.Finalize(skip_zeros);
339 b.Assemble();
340
341 // 7. Define the initial conditions, save the corresponding grid function to
342 // a file and (optionally) save data in the VisIt format and initialize
343 // GLVis visualization.
344 GridFunction u(&fes);
345 u.ProjectCoefficient(u0);
346
347 {
348 ofstream omesh("ex9.mesh");
349 omesh.precision(precision);
350 mesh->Print(omesh);
351 ofstream osol("ex9-init.gf");
352 osol.precision(precision);
353 u.Save(osol);
354 }
355
356 // Create data collection for solution output: either VisItDataCollection for
357 // ascii data files, or SidreDataCollection for binary data files.
358 DataCollection *dc = NULL;
359 if (visit)
360 {
361 if (binary)
362 {
363 #ifdef MFEM_USE_SIDRE
364 dc = new SidreDataCollection("Example9", mesh);
365 #else
366 MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
367 #endif
368 }
369 else
370 {
371 dc = new VisItDataCollection("Example9", mesh);
372 dc->SetPrecision(precision);
373 }
374 dc->RegisterField("solution", &u);
375 dc->SetCycle(0);
376 dc->SetTime(0.0);
377 dc->Save();
378 }
379
380 socketstream sout;
381 if (visualization)
382 {
383 char vishost[] = "localhost";
384 int visport = 19916;
385 sout.open(vishost, visport);
386 if (!sout)
387 {
388 cout << "Unable to connect to GLVis server at "
389 << vishost << ':' << visport << endl;
390 visualization = false;
391 cout << "GLVis visualization disabled.\n";
392 }
393 else
394 {
395 sout.precision(precision);
396 sout << "solution\n" << *mesh << u;
397 sout << "pause\n";
398 sout << flush;
399 cout << "GLVis visualization paused."
400 << " Press space (in the GLVis window) to resume it.\n";
401 }
402 }
403
404 Vector M_rowsums(m.Size());
405 m.SpMat().GetRowSums(M_rowsums);
406
407 // 8. Define the time-dependent evolution operator describing the ODE
408 // right-hand side, and perform time-integration (looping over the time
409 // iterations, ti, with a time-step dt).
410 FE_Evolution adv(m.SpMat(), k.SpMat(), b, k, M_rowsums);
411
412 double t = 0.0;
413 adv.SetTime(t);
414 ode_solver->Init(adv);
415
416 // Compute initial volume.
417 const double vol0 = M_rowsums * u;
418
419 bool done = false;
420 for (int ti = 0; !done; )
421 {
422 double dt_real = min(dt, t_final - t);
423 adv.SetTimeStep(dt_real);
424 ode_solver->Step(u, t, dt_real);
425 ti++;
426
427 done = (t >= t_final - 1e-8*dt);
428
429 if (done || ti % vis_steps == 0)
430 {
431 cout << "time step: " << ti << ", time: " << t << endl;
432
433 if (visualization)
434 {
435 sout << "solution\n" << *mesh << u << flush;
436 }
437
438 if (visit)
439 {
440 dc->SetCycle(ti);
441 dc->SetTime(t);
442 dc->Save();
443 }
444 }
445 }
446
447 // Print the error vs exact solution.
448 const double max_error = u.ComputeMaxError(u0),
449 l1_error = u.ComputeL1Error(u0),
450 l2_error = u.ComputeL2Error(u0);
451 std::cout << "Linf error = " << max_error << endl
452 << "L1 error = " << l1_error << endl
453 << "L2 error = " << l2_error << endl;
454
455 // Print error in volume.
456 const double vol = M_rowsums * u;
457 std::cout << "Vol error = " << vol - vol0 << endl;
458
459 // 9. Save the final solution. This output can be viewed later using GLVis:
460 // "glvis -m ex9.mesh -g ex9-final.gf".
461 {
462 ofstream osol("ex9-final.gf");
463 osol.precision(precision);
464 u.Save(osol);
465 }
466
467 // 10. Free the used memory.
468 delete ode_solver;
469 delete dc;
470 delete mesh;
471
472 return 0;
473 }
474
475
476 // Implementation of class FE_Evolution
FE_Evolution(SparseMatrix & M_,SparseMatrix & K_,const Vector & b_,BilinearForm & bf_,Vector & M_rs)477 FE_Evolution::FE_Evolution(SparseMatrix &M_, SparseMatrix &K_,
478 const Vector &b_, BilinearForm &bf_, Vector &M_rs)
479 : TimeDependentOperator(M_.Size()),
480 M(M_), K(K_), b(b_), M_prec(), M_solver(), z(M_.Size()),
481 bf(bf_), M_rowsums(M_rs)
482 {
483 M_solver.SetPreconditioner(M_prec);
484 M_solver.SetOperator(M);
485
486 M_solver.iterative_mode = false;
487 M_solver.SetRelTol(1e-9);
488 M_solver.SetAbsTol(0.0);
489 M_solver.SetMaxIter(100);
490 M_solver.SetPrintLevel(0);
491 }
492
Mult(const Vector & x,Vector & y) const493 void FE_Evolution::Mult(const Vector &x, Vector &y) const
494 {
495 // Compute bounds y_min, y_max for y from x on the ldofs.
496 const int dofs = x.Size();
497 Vector y_min(dofs), y_max(dofs);
498 const int *In = bf.SpMat().GetI(), *Jn = bf.SpMat().GetJ();
499 for (int i = 0, k = 0; i < dofs; i++)
500 {
501 double x_i_min = +std::numeric_limits<double>::infinity();
502 double x_i_max = -std::numeric_limits<double>::infinity();
503 for (int end = In[i+1]; k < end; k++)
504 {
505 const int j = Jn[k];
506 if (x(j) > x_i_max) { x_i_max = x(j); }
507 if (x(j) < x_i_min) { x_i_min = x(j); }
508 }
509 y_min(i) = x_i_min;
510 y_max(i) = x_i_max;
511 }
512 for (int i = 0; i < dofs; i++)
513 {
514 y_min(i) = (y_min(i) - x(i) ) / dt;
515 y_max(i) = (y_max(i) - x(i) ) / dt;
516 }
517
518 // Compute the high-order solution y = M^{-1} (K x + b).
519 K.Mult(x, z);
520 z += b;
521 M_solver.Mult(z, y);
522
523 // The solution y is an increment; it should not introduce new mass.
524 const double mass_y = 0.0;
525
526 // Perform optimization.
527 Vector y_out(dofs);
528 const int max_iter = 500;
529 const double rtol = 1.e-7;
530 double atol = 1.e-7;
531
532 OptimizationSolver *optsolver = NULL;
533 if (optimizer_type == 2)
534 {
535 #ifdef MFEM_USE_HIOP
536 HiopNlpOptimizer *tmp_opt_ptr = new HiopNlpOptimizer();
537 optsolver = tmp_opt_ptr;
538 #else
539 MFEM_ABORT("MFEM is not built with HiOp support!");
540 #endif
541 }
542 else
543 {
544 SLBQPOptimizer *slbqp = new SLBQPOptimizer();
545 slbqp->SetBounds(y_min, y_max);
546 slbqp->SetLinearConstraint(M_rowsums, mass_y);
547 atol = 1.e-15;
548 optsolver = slbqp;
549 }
550
551 OptimizedTransportProblem ot_prob(y, M_rowsums, mass_y, y_min, y_max);
552 optsolver->SetOptimizationProblem(ot_prob);
553
554 optsolver->SetMaxIter(max_iter);
555 optsolver->SetAbsTol(atol);
556 optsolver->SetRelTol(rtol);
557 optsolver->SetPrintLevel(0);
558 optsolver->Mult(y, y_out);
559
560 y = y_out;
561
562 delete optsolver;
563 }
564
565
566 // Velocity coefficient
velocity_function(const Vector & x,Vector & v)567 void velocity_function(const Vector &x, Vector &v)
568 {
569 int dim = x.Size();
570
571 // map to the reference [-1,1] domain
572 Vector X(dim);
573 for (int i = 0; i < dim; i++)
574 {
575 double center = (bb_min[i] + bb_max[i]) * 0.5;
576 X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
577 }
578
579 switch (problem)
580 {
581 case 0:
582 {
583 // Translations in 1D, 2D, and 3D
584 switch (dim)
585 {
586 case 1: v(0) = (invert_velocity) ? -1.0 : 1.0; break;
587 case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
588 case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
589 break;
590 }
591 break;
592 }
593 case 1:
594 case 2:
595 {
596 // Clockwise rotation in 2D around the origin
597 const double w = M_PI/2;
598 switch (dim)
599 {
600 case 1: v(0) = 1.0; break;
601 case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
602 case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
603 }
604 break;
605 }
606 case 3:
607 {
608 // Clockwise twisting rotation in 2D around the origin
609 const double w = M_PI/2;
610 double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
611 d = d*d;
612 switch (dim)
613 {
614 case 1: v(0) = 1.0; break;
615 case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
616 case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
617 }
618 break;
619 }
620 }
621 }
622
623 // Initial condition
u0_function(const Vector & x)624 double u0_function(const Vector &x)
625 {
626 int dim = x.Size();
627
628 // map to the reference [-1,1] domain
629 Vector X(dim);
630 for (int i = 0; i < dim; i++)
631 {
632 double center = (bb_min[i] + bb_max[i]) * 0.5;
633 X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
634 }
635
636 switch (problem)
637 {
638 case 0:
639 case 1:
640 {
641 switch (dim)
642 {
643 case 1:
644 return (X(0) > -0.15 && X(0) < 0.15) ? 1.0 : 0.0;
645 //return exp(-40.*pow(X(0)-0.0,2));
646 case 2:
647 case 3:
648 {
649 double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
650 if (dim == 3)
651 {
652 const double s = (1. + 0.25*cos(2*M_PI*X(2)));
653 rx *= s;
654 ry *= s;
655 }
656 return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
657 erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
658 }
659 }
660 }
661 case 2:
662 {
663 double x_ = X(0), y_ = X(1), rho, phi;
664 rho = hypot(x_, y_);
665 phi = atan2(y_, x_);
666 return pow(sin(M_PI*rho),2)*sin(3*phi);
667 }
668 case 3:
669 {
670 const double f = M_PI;
671 return sin(f*X(0))*sin(f*X(1));
672 }
673 }
674 return 0.0;
675 }
676
677 // Inflow boundary condition (zero for the problems considered in this example)
inflow_function(const Vector & x)678 double inflow_function(const Vector &x)
679 {
680 switch (problem)
681 {
682 case 0:
683 case 1:
684 case 2:
685 case 3: return 0.0;
686 }
687 return 0.0;
688 }
689