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