MFEM  v3.4 Finite element discretization library
ex9p.cpp
Go to the documentation of this file.
1 // MFEM Example 9 - Parallel Version
2 // SUNDIALS Modification
3 //
4 // Compile with: make ex9p
5 //
6 // Sample runs:
7 // mpirun -np 4 ex9p -m ../../data/periodic-segment.mesh -p 1 -rp 1 -s 11 -dt 0.0025
8 // mpirun -np 4 ex9p -m ../../data/periodic-square.mesh -p 1 -rp 1 -s 12 -dt 0.0025 -tf 9
9 // mpirun -np 4 ex9p -m ../../data/periodic-hexagon.mesh -p 0 -rp 1 -s 11 -dt 0.0009 -vs 25
10 // mpirun -np 4 ex9p -m ../../data/periodic-hexagon.mesh -p 0 -rp 1 -s 13 -dt 0.005 -vs 15
11 // mpirun -np 4 ex9p -m ../../data/amr-quad.mesh -p 1 -rp 1 -s 13 -dt 0.001 -tf 9
12 // mpirun -np 4 ex9p -m ../../data/star-q3.mesh -p 1 -rp 1 -s 13 -dt 0.0025 -tf 9
13 // mpirun -np 4 ex9p -m ../../data/disc-nurbs.mesh -p 1 -rp 2 -s 11 -dt 0.0025 -tf 9
14 // mpirun -np 4 ex9p -m ../../data/periodic-cube.mesh -p 0 -rp 1 -s 12 -dt 0.01 -tf 8 -o 2
15 //
16 // Description: This example code solves the time-dependent advection equation
17 // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
18 // u0(x)=u(0,x) is a given initial condition.
19 //
20 // The example demonstrates the use of Discontinuous Galerkin (DG)
21 // bilinear forms in MFEM (face integrators), the use of explicit
22 // ODE time integrators, the definition of periodic boundary
23 // conditions through periodic meshes, as well as the use of GLVis
24 // for persistent visualization of a time-evolving solution. The
25 // saving of time-dependent data files for external visualization
26 // with VisIt (visit.llnl.gov) is also illustrated.
27
28 #include "mfem.hpp"
29 #include <fstream>
30 #include <iostream>
31
32 #ifndef MFEM_USE_SUNDIALS
33 #error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
34 #endif
35
36 using namespace std;
37 using namespace mfem;
38
39 // Choice for the problem setup. The fluid velocity, initial condition and
40 // inflow boundary condition are chosen based on this parameter.
41 int problem;
42
43 // Velocity coefficient
44 void velocity_function(const Vector &x, Vector &v);
45
46 // Initial condition
47 double u0_function(const Vector &x);
48
49 // Inflow boundary condition
50 double inflow_function(const Vector &x);
51
52 // Mesh bounding box
54
55
56 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
57  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
58  and advection matrices, and b describes the flow on the boundary. This can
59  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
60  used to evaluate the right-hand side. */
62 {
63 private:
64  HypreParMatrix &M, &K;
65  const Vector &b;
66  HypreSmoother M_prec;
67  CGSolver M_solver;
68
69  mutable Vector z;
70
71 public:
72  FE_Evolution(HypreParMatrix &_M, HypreParMatrix &_K, const Vector &_b);
73
74  virtual void Mult(const Vector &x, Vector &y) const;
75
76  virtual ~FE_Evolution() { }
77 };
78
79
80 int main(int argc, char *argv[])
81 {
82  // 1. Initialize MPI.
83  int num_procs, myid;
84  MPI_Init(&argc, &argv);
85  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
86  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
87
88  // 2. Parse command-line options.
89  problem = 0;
90  const char *mesh_file = "../../data/periodic-hexagon.mesh";
91  int ser_ref_levels = 2;
92  int par_ref_levels = 0;
93  int order = 3;
94  int ode_solver_type = 4;
95  double t_final = 10.0;
96  double dt = 0.01;
97  bool visualization = true;
98  bool visit = false;
99  bool binary = false;
100  int vis_steps = 5;
101
102  // Relative and absolute tolerances for CVODE and ARKODE.
103  const double reltol = 1e-2, abstol = 1e-2;
104
105  int precision = 8;
106  cout.precision(precision);
107
108  OptionsParser args(argc, argv);
110  "Mesh file to use.");
112  "Problem setup to use. See options in velocity_function().");
114  "Number of times to refine the mesh uniformly in serial.");
116  "Number of times to refine the mesh uniformly in parallel.");
118  "Order (degree) of the finite elements.");
120  "ODE solver: 1 - Forward Euler,\n\t"
121  " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6,\n\t"
122  " 11 - CVODE (adaptive order) explicit,\n\t"
123  " 12 - ARKODE default (4th order) explicit,\n\t"
124  " 13 - ARKODE RK8.");
126  "Final time; start time is 0.");
128  "Time step.");
130  "--no-visualization",
131  "Enable or disable GLVis visualization.");
133  "--no-visit-datafiles",
134  "Save data files for VisIt (visit.llnl.gov) visualization.");
136  "--ascii-datafiles",
137  "Use binary (Sidre) or ascii format for VisIt data files.");
139  "Visualize every n-th timestep.");
140  args.Parse();
141  if (!args.Good())
142  {
143  if (myid == 0)
144  {
145  args.PrintUsage(cout);
146  }
147  MPI_Finalize();
148  return 1;
149  }
150  if (myid == 0)
151  {
152  args.PrintOptions(cout);
153  }
154
155  // 3. Read the serial mesh from the given mesh file on all processors. We can
156  // handle geometrically periodic meshes in this code.
157  Mesh *mesh = new Mesh(mesh_file, 1, 1);
158  int dim = mesh->Dimension();
159
160  // 4. Define the ODE solver used for time integration. Several explicit
161  // Runge-Kutta methods are available.
162  ODESolver *ode_solver = NULL;
163  CVODESolver *cvode = NULL;
164  ARKODESolver *arkode = NULL;
165  switch (ode_solver_type)
166  {
167  case 1: ode_solver = new ForwardEulerSolver; break;
168  case 2: ode_solver = new RK2Solver(1.0); break;
169  case 3: ode_solver = new RK3SSPSolver; break;
170  case 4: ode_solver = new RK4Solver; break;
171  case 6: ode_solver = new RK6Solver; break;
172  case 11:
173  cvode = new CVODESolver(MPI_COMM_WORLD, CV_ADAMS, CV_FUNCTIONAL);
174  cvode->SetSStolerances(reltol, abstol);
175  cvode->SetMaxStep(dt);
176  ode_solver = cvode; break;
177  case 12:
178  case 13:
179  arkode = new ARKODESolver(MPI_COMM_WORLD, ARKODESolver::EXPLICIT);
180  arkode->SetSStolerances(reltol, abstol);
181  arkode->SetMaxStep(dt);
182  if (ode_solver_type == 13) { arkode->SetERKTableNum(FEHLBERG_13_7_8); }
183  ode_solver = arkode; break;
184  default:
185  if (myid == 0)
186  {
187  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
188  }
189  delete mesh;
190  MPI_Finalize();
191  return 3;
192  }
193
194  // 5. Refine the mesh in serial to increase the resolution. In this example
195  // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
196  // a command-line parameter. If the mesh is of NURBS type, we convert it
197  // to a (piecewise-polynomial) high-order mesh.
198  for (int lev = 0; lev < ser_ref_levels; lev++)
199  {
200  mesh->UniformRefinement();
201  }
202  if (mesh->NURBSext)
203  {
204  mesh->SetCurvature(max(order, 1));
205  }
206  mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
207
208  // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
209  // this mesh further in parallel to increase the resolution. Once the
210  // parallel mesh is defined, the serial mesh can be deleted.
211  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
212  delete mesh;
213  for (int lev = 0; lev < par_ref_levels; lev++)
214  {
215  pmesh->UniformRefinement();
216  }
217
218  // 7. Define the parallel discontinuous DG finite element space on the
219  // parallel refined mesh of the given polynomial order.
220  DG_FECollection fec(order, dim);
221  ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
222
223  HYPRE_Int global_vSize = fes->GlobalTrueVSize();
224  if (myid == 0)
225  {
226  cout << "Number of unknowns: " << global_vSize << endl;
227  }
228
229  // 8. Set up and assemble the parallel bilinear and linear forms (and the
230  // parallel hypre matrices) corresponding to the DG discretization. The
231  // DGTraceIntegrator involves integrals over mesh interior faces.
235
236  ParBilinearForm *m = new ParBilinearForm(fes);
238  ParBilinearForm *k = new ParBilinearForm(fes);
241  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
243  new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
244
245  ParLinearForm *b = new ParLinearForm(fes);
247  new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
248
249  m->Assemble();
250  m->Finalize();
251  int skip_zeros = 0;
252  k->Assemble(skip_zeros);
253  k->Finalize(skip_zeros);
254  b->Assemble();
255
259
260  // 9. Define the initial conditions, save the corresponding grid function to
261  // a file and (optionally) save data in the VisIt format and initialize
262  // GLVis visualization.
263  ParGridFunction *u = new ParGridFunction(fes);
264  u->ProjectCoefficient(u0);
265  HypreParVector *U = u->GetTrueDofs();
266
267  {
268  ostringstream mesh_name, sol_name;
269  mesh_name << "ex9-mesh." << setfill('0') << setw(6) << myid;
270  sol_name << "ex9-init." << setfill('0') << setw(6) << myid;
271  ofstream omesh(mesh_name.str().c_str());
272  omesh.precision(precision);
273  pmesh->Print(omesh);
274  ofstream osol(sol_name.str().c_str());
275  osol.precision(precision);
276  u->Save(osol);
277  }
278
279  // Create data collection for solution output: either VisItDataCollection for
280  // ascii data files, or SidreDataCollection for binary data files.
281  DataCollection *dc = NULL;
282  if (visit)
283  {
284  if (binary)
285  {
286 #ifdef MFEM_USE_SIDRE
287  dc = new SidreDataCollection("Example9-Parallel", pmesh);
288 #else
289  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
290 #endif
291  }
292  else
293  {
294  dc = new VisItDataCollection("Example9-Parallel", pmesh);
295  dc->SetPrecision(precision);
296  }
297  dc->RegisterField("solution", u);
298  dc->SetCycle(0);
299  dc->SetTime(0.0);
300  dc->Save();
301  }
302
303  socketstream sout;
304  if (visualization)
305  {
306  char vishost[] = "localhost";
307  int visport = 19916;
308  sout.open(vishost, visport);
309  if (!sout)
310  {
311  if (myid == 0)
312  cout << "Unable to connect to GLVis server at "
313  << vishost << ':' << visport << endl;
314  visualization = false;
315  if (myid == 0)
316  {
317  cout << "GLVis visualization disabled.\n";
318  }
319  }
320  else
321  {
322  sout << "parallel " << num_procs << " " << myid << "\n";
323  sout.precision(precision);
324  sout << "solution\n" << *pmesh << *u;
325  sout << "pause\n";
326  sout << flush;
327  if (myid == 0)
328  cout << "GLVis visualization paused."
329  << " Press space (in the GLVis window) to resume it.\n";
330  }
331  }
332
333  // 10. Define the time-dependent evolution operator describing the ODE
334  // right-hand side, and perform time-integration (looping over the time
335  // iterations, ti, with a time-step dt).
337
338  double t = 0.0;
341
342  bool done = false;
343  for (int ti = 0; !done; )
344  {
345  double dt_real = min(dt, t_final - t);
346  ode_solver->Step(*U, t, dt_real);
347  ti++;
348
349  done = (t >= t_final - 1e-8*dt);
350
351  if (done || ti % vis_steps == 0)
352  {
353  if (myid == 0)
354  {
355  cout << "time step: " << ti << ", time: " << t << endl;
356  if (cvode) { cvode->PrintInfo(); }
357  if (arkode) { arkode->PrintInfo(); }
358  }
359
360  // 11. Extract the parallel grid function corresponding to the finite
361  // element approximation U (the local solution on each processor).
362  *u = *U;
363
364  if (visualization)
365  {
366  sout << "parallel " << num_procs << " " << myid << "\n";
367  sout << "solution\n" << *pmesh << *u << flush;
368  }
369
370  if (visit)
371  {
372  dc->SetCycle(ti);
373  dc->SetTime(t);
374  dc->Save();
375  }
376  }
377  }
378
379  // 12. Save the final solution in parallel. This output can be viewed later
380  // using GLVis: "glvis -np <np> -m ex9-mesh -g ex9-final".
381  {
382  *u = *U;
383  ostringstream sol_name;
384  sol_name << "ex9-final." << setfill('0') << setw(6) << myid;
385  ofstream osol(sol_name.str().c_str());
386  osol.precision(precision);
387  u->Save(osol);
388  }
389
390  // 13. Free the used memory.
391  delete U;
392  delete u;
393  delete B;
394  delete b;
395  delete K;
396  delete k;
397  delete M;
398  delete m;
399  delete fes;
400  delete pmesh;
401  delete ode_solver;
402  delete dc;
403
404  MPI_Finalize();
405  return 0;
406 }
407
408
409 // Implementation of class FE_Evolution
411  const Vector &_b)
412  : TimeDependentOperator(_M.Height()),
413  M(_M), K(_K), b(_b), M_solver(M.GetComm()), z(_M.Height())
414 {
415  M_prec.SetType(HypreSmoother::Jacobi);
416  M_solver.SetPreconditioner(M_prec);
417  M_solver.SetOperator(M);
418
419  M_solver.iterative_mode = false;
420  M_solver.SetRelTol(1e-9);
421  M_solver.SetAbsTol(0.0);
422  M_solver.SetMaxIter(100);
423  M_solver.SetPrintLevel(0);
424 }
425
426 void FE_Evolution::Mult(const Vector &x, Vector &y) const
427 {
428  // y = M^{-1} (K x + b)
429  K.Mult(x, z);
430  z += b;
431  M_solver.Mult(z, y);
432 }
433
434
435 // Velocity coefficient
436 void velocity_function(const Vector &x, Vector &v)
437 {
438  int dim = x.Size();
439
440  // map to the reference [-1,1] domain
441  Vector X(dim);
442  for (int i = 0; i < dim; i++)
443  {
444  double center = (bb_min[i] + bb_max[i]) * 0.5;
445  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
446  }
447
448  switch (problem)
449  {
450  case 0:
451  {
452  // Translations in 1D, 2D, and 3D
453  switch (dim)
454  {
455  case 1: v(0) = 1.0; break;
456  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
457  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
458  break;
459  }
460  break;
461  }
462  case 1:
463  case 2:
464  {
465  // Clockwise rotation in 2D around the origin
466  const double w = M_PI/2;
467  switch (dim)
468  {
469  case 1: v(0) = 1.0; break;
470  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
471  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
472  }
473  break;
474  }
475  case 3:
476  {
477  // Clockwise twisting rotation in 2D around the origin
478  const double w = M_PI/2;
479  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
480  d = d*d;
481  switch (dim)
482  {
483  case 1: v(0) = 1.0; break;
484  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
485  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
486  }
487  break;
488  }
489  }
490 }
491
492 // Initial condition
493 double u0_function(const Vector &x)
494 {
495  int dim = x.Size();
496
497  // map to the reference [-1,1] domain
498  Vector X(dim);
499  for (int i = 0; i < dim; i++)
500  {
501  double center = (bb_min[i] + bb_max[i]) * 0.5;
502  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
503  }
504
505  switch (problem)
506  {
507  case 0:
508  case 1:
509  {
510  switch (dim)
511  {
512  case 1:
513  return exp(-40.*pow(X(0)-0.5,2));
514  case 2:
515  case 3:
516  {
517  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
518  if (dim == 3)
519  {
520  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
521  rx *= s;
522  ry *= s;
523  }
524  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
525  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
526  }
527  }
528  }
529  case 2:
530  {
531  double x_ = X(0), y_ = X(1), rho, phi;
532  rho = hypot(x_, y_);
533  phi = atan2(y_, x_);
534  return pow(sin(M_PI*rho),2)*sin(3*phi);
535  }
536  case 3:
537  {
538  const double f = M_PI;
539  return sin(f*X(0))*sin(f*X(1));
540  }
541  }
542  return 0.0;
543 }
544
545 // Inflow boundary condition (zero for the problems considered in this example)
546 double inflow_function(const Vector &x)
547 {
548  switch (problem)
549  {
550  case 0:
551  case 1:
552  case 2:
553  case 3: return 0.0;
554  }
555  return 0.0;
556 }
Vector bb_max
Definition: ex9.cpp:52
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
Definition: solvers.hpp:111
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:224
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:295
virtual void Init(TimeDependentOperator &f)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.hpp:37
Base abstract class for time dependent operators.
Definition: operator.hpp:151
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Definition: mesh.cpp:108
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:51
virtual void Step(Vector &x, double &t, double &dt)=0
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in]...
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:181
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:398
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:290
int main(int argc, char *argv[])
Definition: ex1.cpp:45
void SetSStolerances(double reltol, double abstol)
Specify the scalar relative and scalar absolute tolerances.
Definition: sundials.cpp:477
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void SetERKTableNum(int table_num)
Choose a specific Butcher table for explicit RK method.
Definition: sundials.cpp:528
FE_Evolution(FiniteElementSpace &_vfes, Operator &_A, SparseMatrix &_Aflux)
Definition: ex18.hpp:104
Data collection with Sidre routines following the Conduit mesh blueprint specification.
Class for parallel linear form.
Definition: plinearform.hpp:26
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
int dim
Definition: ex3.cpp:47
virtual void Save()
Save the collection to disk.
Wrapper for SUNDIALS&#39; CVODE library – Multi-step time integration.
Definition: sundials.hpp:133
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6741
Data collection with VisIt I/O routines.
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:3355
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:247
void SetMaxStep(double dt_max)
Set the maximum time step of the linear multistep method.
Definition: sundials.hpp:183
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:3428
Parallel smoothers in hypre.
Definition: hypre.hpp:517
Definition: linearform.cpp:38
int Dimension() const
Definition: mesh.hpp:645
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
void SetTime(double t)
Set physical time (for time-dependent simulations)
Wrapper for SUNDIALS&#39; ARKODE library – Runge-Kutta time integration.
Definition: sundials.hpp:220
Vector bb_min
Definition: ex9.cpp:52
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:73
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:147
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:340
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:486
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Definition: optparser.hpp:74
int problem
Definition: ex15.cpp:57
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:134
virtual ~FE_Evolution()
Definition: ex9p.cpp:76
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: ex18.hpp:130
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:176
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:141
Class for parallel bilinear form.
int open(const char hostname[], int port)
double u0_function(const Vector &x)
Definition: ex9.cpp:397
class for C-function coefficient
Vector data type.
Definition: vector.hpp:48
void PrintInfo() const
Print CVODE statistics.
Definition: sundials.cpp:400
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The classical forward Euler method.
Definition: ode.hpp:101
void SetMaxStep(double dt_max)
Set the maximum time step of the Runge-Kutta method.
Definition: sundials.hpp:276
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
Class for parallel meshes.
Definition: pmesh.hpp:32
void PrintInfo() const
Print ARKODE statistics.
Definition: sundials.cpp:691
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:34
double inflow_function(const Vector &x)
Definition: ex9.cpp:450