MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex9.cpp
Go to the documentation of this file.
1 // MFEM Example 9
2 //
3 // Compile with: make ex9
4 //
5 // Sample runs:
6 // ex9 -m ../data/periodic-segment.mesh -p 0 -r 2 -dt 0.005
7 // ex9 -m ../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10
8 // ex9 -m ../data/periodic-hexagon.mesh -p 0 -r 2 -dt 0.01 -tf 10
9 // ex9 -m ../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9
10 // ex9 -m ../data/periodic-hexagon.mesh -p 1 -r 2 -dt 0.005 -tf 9
11 // ex9 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9
12 // ex9 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.02 -s 13 -tf 9
13 // ex9 -m ../data/star-q3.mesh -p 1 -r 2 -dt 0.005 -tf 9
14 // ex9 -m ../data/star-mixed.mesh -p 1 -r 2 -dt 0.005 -tf 9
15 // ex9 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.005 -tf 9
16 // ex9 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9
17 // ex9 -m ../data/periodic-square.mesh -p 3 -r 4 -dt 0.0025 -tf 9 -vs 20
18 // ex9 -m ../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8
19 // ex9 -m ../data/periodic-square.msh -p 0 -r 2 -dt 0.005 -tf 2
20 // ex9 -m ../data/periodic-cube.msh -p 0 -r 1 -o 2 -tf 2
21 //
22 // Device sample runs:
23 // ex9 -pa
24 // ex9 -ea
25 // ex9 -fa
26 // ex9 -pa -m ../data/periodic-cube.mesh
27 // ex9 -pa -m ../data/periodic-cube.mesh -d cuda
28 // ex9 -ea -m ../data/periodic-cube.mesh -d cuda
29 // ex9 -fa -m ../data/periodic-cube.mesh -d cuda
30 // ex9 -pa -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9 -d cuda
31 //
32 // Description: This example code solves the time-dependent advection equation
33 // du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
34 // u0(x)=u(0,x) is a given initial condition.
35 //
36 // The example demonstrates the use of Discontinuous Galerkin (DG)
37 // bilinear forms in MFEM (face integrators), the use of implicit
38 // and explicit ODE time integrators, the definition of periodic
39 // boundary conditions through periodic meshes, as well as the use
40 // of GLVis for persistent visualization of a time-evolving
41 // solution. The saving of time-dependent data files for external
42 // visualization with VisIt (visit.llnl.gov) and ParaView
43 // (paraview.org) is also illustrated.
44 
45 #include "mfem.hpp"
46 #include <fstream>
47 #include <iostream>
48 #include <algorithm>
49 
50 using namespace std;
51 using namespace mfem;
52 
53 // Choice for the problem setup. The fluid velocity, initial condition and
54 // inflow boundary condition are chosen based on this parameter.
55 int problem;
56 
57 // Velocity coefficient
58 void velocity_function(const Vector &x, Vector &v);
59 
60 // Initial condition
61 double u0_function(const Vector &x);
62 
63 // Inflow boundary condition
64 double inflow_function(const Vector &x);
65 
66 // Mesh bounding box
68 
69 class DG_Solver : public Solver
70 {
71 private:
72  SparseMatrix &M, &K, A;
73  GMRESSolver linear_solver;
74  BlockILU prec;
75  double dt;
76 public:
77  DG_Solver(SparseMatrix &M_, SparseMatrix &K_, const FiniteElementSpace &fes)
78  : M(M_),
79  K(K_),
80  prec(fes.GetFE(0)->GetDof(),
81  BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
82  dt(-1.0)
83  {
84  linear_solver.iterative_mode = false;
85  linear_solver.SetRelTol(1e-9);
86  linear_solver.SetAbsTol(0.0);
87  linear_solver.SetMaxIter(100);
88  linear_solver.SetPrintLevel(0);
89  linear_solver.SetPreconditioner(prec);
90  }
91 
92  void SetTimeStep(double dt_)
93  {
94  if (dt_ != dt)
95  {
96  dt = dt_;
97  // Form operator A = M - dt*K
98  A = K;
99  A *= -dt;
100  A += M;
101 
102  // this will also call SetOperator on the preconditioner
103  linear_solver.SetOperator(A);
104  }
105  }
106 
107  void SetOperator(const Operator &op)
108  {
109  linear_solver.SetOperator(op);
110  }
111 
112  virtual void Mult(const Vector &x, Vector &y) const
113  {
114  linear_solver.Mult(x, y);
115  }
116 };
117 
118 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
119  form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
120  and advection matrices, and b describes the flow on the boundary. This can
121  be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
122  used to evaluate the right-hand side. */
123 class FE_Evolution : public TimeDependentOperator
124 {
125 private:
126  BilinearForm &M, &K;
127  const Vector &b;
128  Solver *M_prec;
129  CGSolver M_solver;
130  DG_Solver *dg_solver;
131 
132  mutable Vector z;
133 
134 public:
135  FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_);
136 
137  virtual void Mult(const Vector &x, Vector &y) const;
138  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
139 
140  virtual ~FE_Evolution();
141 };
142 
143 
144 int main(int argc, char *argv[])
145 {
146  // 1. Parse command-line options.
147  problem = 0;
148  const char *mesh_file = "../data/periodic-hexagon.mesh";
149  int ref_levels = 2;
150  int order = 3;
151  bool pa = false;
152  bool ea = false;
153  bool fa = false;
154  const char *device_config = "cpu";
155  int ode_solver_type = 4;
156  double t_final = 10.0;
157  double dt = 0.01;
158  bool visualization = true;
159  bool visit = false;
160  bool paraview = false;
161  bool binary = false;
162  int vis_steps = 5;
163 
164  int precision = 8;
165  cout.precision(precision);
166 
167  OptionsParser args(argc, argv);
168  args.AddOption(&mesh_file, "-m", "--mesh",
169  "Mesh file to use.");
170  args.AddOption(&problem, "-p", "--problem",
171  "Problem setup to use. See options in velocity_function().");
172  args.AddOption(&ref_levels, "-r", "--refine",
173  "Number of times to refine the mesh uniformly.");
174  args.AddOption(&order, "-o", "--order",
175  "Order (degree) of the finite elements.");
176  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
177  "--no-partial-assembly", "Enable Partial Assembly.");
178  args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
179  "--no-element-assembly", "Enable Element Assembly.");
180  args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
181  "--no-full-assembly", "Enable Full Assembly.");
182  args.AddOption(&device_config, "-d", "--device",
183  "Device configuration string, see Device::Configure().");
184  args.AddOption(&ode_solver_type, "-s", "--ode-solver",
185  "ODE solver: 1 - Forward Euler,\n\t"
186  " 2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6,\n\t"
187  " 11 - Backward Euler,\n\t"
188  " 12 - SDIRK23 (L-stable), 13 - SDIRK33,\n\t"
189  " 22 - Implicit Midpoint Method,\n\t"
190  " 23 - SDIRK23 (A-stable), 24 - SDIRK34");
191  args.AddOption(&t_final, "-tf", "--t-final",
192  "Final time; start time is 0.");
193  args.AddOption(&dt, "-dt", "--time-step",
194  "Time step.");
195  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
196  "--no-visualization",
197  "Enable or disable GLVis visualization.");
198  args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
199  "--no-visit-datafiles",
200  "Save data files for VisIt (visit.llnl.gov) visualization.");
201  args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
202  "--no-paraview-datafiles",
203  "Save data files for ParaView (paraview.org) visualization.");
204  args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
205  "--ascii-datafiles",
206  "Use binary (Sidre) or ascii format for VisIt data files.");
207  args.AddOption(&vis_steps, "-vs", "--visualization-steps",
208  "Visualize every n-th timestep.");
209  args.Parse();
210  if (!args.Good())
211  {
212  args.PrintUsage(cout);
213  return 1;
214  }
215  args.PrintOptions(cout);
216 
217  Device device(device_config);
218  device.Print();
219 
220  // 2. Read the mesh from the given mesh file. We can handle geometrically
221  // periodic meshes in this code.
222  Mesh mesh(mesh_file, 1, 1);
223  int dim = mesh.Dimension();
224 
225  // 3. Define the ODE solver used for time integration. Several explicit
226  // Runge-Kutta methods are available.
227  ODESolver *ode_solver = NULL;
228  switch (ode_solver_type)
229  {
230  // Explicit methods
231  case 1: ode_solver = new ForwardEulerSolver; break;
232  case 2: ode_solver = new RK2Solver(1.0); break;
233  case 3: ode_solver = new RK3SSPSolver; break;
234  case 4: ode_solver = new RK4Solver; break;
235  case 6: ode_solver = new RK6Solver; break;
236  // Implicit (L-stable) methods
237  case 11: ode_solver = new BackwardEulerSolver; break;
238  case 12: ode_solver = new SDIRK23Solver(2); break;
239  case 13: ode_solver = new SDIRK33Solver; break;
240  // Implicit A-stable methods (not L-stable)
241  case 22: ode_solver = new ImplicitMidpointSolver; break;
242  case 23: ode_solver = new SDIRK23Solver; break;
243  case 24: ode_solver = new SDIRK34Solver; break;
244 
245  default:
246  cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
247  return 3;
248  }
249 
250  // 4. Refine the mesh to increase the resolution. In this example we do
251  // 'ref_levels' of uniform refinement, where 'ref_levels' is a
252  // command-line parameter. If the mesh is of NURBS type, we convert it to
253  // a (piecewise-polynomial) high-order mesh.
254  for (int lev = 0; lev < ref_levels; lev++)
255  {
256  mesh.UniformRefinement();
257  }
258  if (mesh.NURBSext)
259  {
260  mesh.SetCurvature(max(order, 1));
261  }
262  mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
263 
264  // 5. Define the discontinuous DG finite element space of the given
265  // polynomial order on the refined mesh.
266  DG_FECollection fec(order, dim, BasisType::GaussLobatto);
267  FiniteElementSpace fes(&mesh, &fec);
268 
269  cout << "Number of unknowns: " << fes.GetVSize() << endl;
270 
271  // 6. Set up and assemble the bilinear and linear forms corresponding to the
272  // DG discretization. The DGTraceIntegrator involves integrals over mesh
273  // interior faces.
277 
278  BilinearForm m(&fes);
279  BilinearForm k(&fes);
280  if (pa)
281  {
282  m.SetAssemblyLevel(AssemblyLevel::PARTIAL);
283  k.SetAssemblyLevel(AssemblyLevel::PARTIAL);
284  }
285  else if (ea)
286  {
287  m.SetAssemblyLevel(AssemblyLevel::ELEMENT);
288  k.SetAssemblyLevel(AssemblyLevel::ELEMENT);
289  }
290  else if (fa)
291  {
292  m.SetAssemblyLevel(AssemblyLevel::FULL);
293  k.SetAssemblyLevel(AssemblyLevel::FULL);
294  }
296  constexpr double alpha = -1.0;
297  k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
299  new NonconservativeDGTraceIntegrator(velocity, alpha));
301  new NonconservativeDGTraceIntegrator(velocity, alpha));
302 
303  LinearForm b(&fes);
305  new BoundaryFlowIntegrator(inflow, velocity, alpha));
306 
307  m.Assemble();
308  int skip_zeros = 0;
309  k.Assemble(skip_zeros);
310  b.Assemble();
311  m.Finalize();
312  k.Finalize(skip_zeros);
313 
314  // 7. Define the initial conditions, save the corresponding grid function to
315  // a file and (optionally) save data in the VisIt format and initialize
316  // GLVis visualization.
317  GridFunction u(&fes);
318  u.ProjectCoefficient(u0);
319 
320  {
321  ofstream omesh("ex9.mesh");
322  omesh.precision(precision);
323  mesh.Print(omesh);
324  ofstream osol("ex9-init.gf");
325  osol.precision(precision);
326  u.Save(osol);
327  }
328 
329  // Create data collection for solution output: either VisItDataCollection for
330  // ascii data files, or SidreDataCollection for binary data files.
331  DataCollection *dc = NULL;
332  if (visit)
333  {
334  if (binary)
335  {
336 #ifdef MFEM_USE_SIDRE
337  dc = new SidreDataCollection("Example9", &mesh);
338 #else
339  MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
340 #endif
341  }
342  else
343  {
344  dc = new VisItDataCollection("Example9", &mesh);
345  dc->SetPrecision(precision);
346  }
347  dc->RegisterField("solution", &u);
348  dc->SetCycle(0);
349  dc->SetTime(0.0);
350  dc->Save();
351  }
352 
353  ParaViewDataCollection *pd = NULL;
354  if (paraview)
355  {
356  pd = new ParaViewDataCollection("Example9", &mesh);
357  pd->SetPrefixPath("ParaView");
358  pd->RegisterField("solution", &u);
359  pd->SetLevelsOfDetail(order);
360  pd->SetDataFormat(VTKFormat::BINARY);
361  pd->SetHighOrderOutput(true);
362  pd->SetCycle(0);
363  pd->SetTime(0.0);
364  pd->Save();
365  }
366 
367  socketstream sout;
368  if (visualization)
369  {
370  char vishost[] = "localhost";
371  int visport = 19916;
372  sout.open(vishost, visport);
373  if (!sout)
374  {
375  cout << "Unable to connect to GLVis server at "
376  << vishost << ':' << visport << endl;
377  visualization = false;
378  cout << "GLVis visualization disabled.\n";
379  }
380  else
381  {
382  sout.precision(precision);
383  sout << "solution\n" << mesh << u;
384  sout << "pause\n";
385  sout << flush;
386  cout << "GLVis visualization paused."
387  << " Press space (in the GLVis window) to resume it.\n";
388  }
389  }
390 
391  // 8. Define the time-dependent evolution operator describing the ODE
392  // right-hand side, and perform time-integration (looping over the time
393  // iterations, ti, with a time-step dt).
394  FE_Evolution adv(m, k, b);
395 
396  double t = 0.0;
397  adv.SetTime(t);
398  ode_solver->Init(adv);
399 
400  bool done = false;
401  for (int ti = 0; !done; )
402  {
403  double dt_real = min(dt, t_final - t);
404  ode_solver->Step(u, t, dt_real);
405  ti++;
406 
407  done = (t >= t_final - 1e-8*dt);
408 
409  if (done || ti % vis_steps == 0)
410  {
411  cout << "time step: " << ti << ", time: " << t << endl;
412 
413  if (visualization)
414  {
415  sout << "solution\n" << mesh << u << flush;
416  }
417 
418  if (visit)
419  {
420  dc->SetCycle(ti);
421  dc->SetTime(t);
422  dc->Save();
423  }
424 
425  if (paraview)
426  {
427  pd->SetCycle(ti);
428  pd->SetTime(t);
429  pd->Save();
430  }
431  }
432  }
433 
434  // 9. Save the final solution. This output can be viewed later using GLVis:
435  // "glvis -m ex9.mesh -g ex9-final.gf".
436  {
437  ofstream osol("ex9-final.gf");
438  osol.precision(precision);
439  u.Save(osol);
440  }
441 
442  // 10. Free the used memory.
443  delete ode_solver;
444  delete pd;
445  delete dc;
446 
447  return 0;
448 }
449 
450 
451 // Implementation of class FE_Evolution
453  : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), z(M_.Height())
454 {
455  Array<int> ess_tdof_list;
456  if (M.GetAssemblyLevel() == AssemblyLevel::LEGACY)
457  {
458  M_prec = new DSmoother(M.SpMat());
459  M_solver.SetOperator(M.SpMat());
460  dg_solver = new DG_Solver(M.SpMat(), K.SpMat(), *M.FESpace());
461  }
462  else
463  {
464  M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
465  M_solver.SetOperator(M);
466  dg_solver = NULL;
467  }
468  M_solver.SetPreconditioner(*M_prec);
469  M_solver.iterative_mode = false;
470  M_solver.SetRelTol(1e-9);
471  M_solver.SetAbsTol(0.0);
472  M_solver.SetMaxIter(100);
473  M_solver.SetPrintLevel(0);
474 }
475 
476 void FE_Evolution::Mult(const Vector &x, Vector &y) const
477 {
478  // y = M^{-1} (K x + b)
479  K.Mult(x, z);
480  z += b;
481  M_solver.Mult(z, y);
482 }
483 
484 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
485 {
486  MFEM_VERIFY(dg_solver != NULL,
487  "Implicit time integration is not supported with partial assembly");
488  K.Mult(x, z);
489  z += b;
490  dg_solver->SetTimeStep(dt);
491  dg_solver->Mult(z, k);
492 }
493 
495 {
496  delete M_prec;
497  delete dg_solver;
498 }
499 
500 // Velocity coefficient
501 void velocity_function(const Vector &x, Vector &v)
502 {
503  int dim = x.Size();
504 
505  // map to the reference [-1,1] domain
506  Vector X(dim);
507  for (int i = 0; i < dim; i++)
508  {
509  double center = (bb_min[i] + bb_max[i]) * 0.5;
510  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
511  }
512 
513  switch (problem)
514  {
515  case 0:
516  {
517  // Translations in 1D, 2D, and 3D
518  switch (dim)
519  {
520  case 1: v(0) = 1.0; break;
521  case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
522  case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
523  break;
524  }
525  break;
526  }
527  case 1:
528  case 2:
529  {
530  // Clockwise rotation in 2D around the origin
531  const double w = M_PI/2;
532  switch (dim)
533  {
534  case 1: v(0) = 1.0; break;
535  case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
536  case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
537  }
538  break;
539  }
540  case 3:
541  {
542  // Clockwise twisting rotation in 2D around the origin
543  const double w = M_PI/2;
544  double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
545  d = d*d;
546  switch (dim)
547  {
548  case 1: v(0) = 1.0; break;
549  case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
550  case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
551  }
552  break;
553  }
554  }
555 }
556 
557 // Initial condition
558 double u0_function(const Vector &x)
559 {
560  int dim = x.Size();
561 
562  // map to the reference [-1,1] domain
563  Vector X(dim);
564  for (int i = 0; i < dim; i++)
565  {
566  double center = (bb_min[i] + bb_max[i]) * 0.5;
567  X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
568  }
569 
570  switch (problem)
571  {
572  case 0:
573  case 1:
574  {
575  switch (dim)
576  {
577  case 1:
578  return exp(-40.*pow(X(0)-0.5,2));
579  case 2:
580  case 3:
581  {
582  double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
583  if (dim == 3)
584  {
585  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
586  rx *= s;
587  ry *= s;
588  }
589  return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
590  erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
591  }
592  }
593  }
594  case 2:
595  {
596  double x_ = X(0), y_ = X(1), rho, phi;
597  rho = hypot(x_, y_);
598  phi = atan2(y_, x_);
599  return pow(sin(M_PI*rho),2)*sin(3*phi);
600  }
601  case 3:
602  {
603  const double f = M_PI;
604  return sin(f*X(0))*sin(f*X(1));
605  }
606  }
607  return 0.0;
608 }
609 
610 // Inflow boundary condition (zero for the problems considered in this example)
611 double inflow_function(const Vector &x)
612 {
613  switch (problem)
614  {
615  case 0:
616  case 1:
617  case 2:
618  case 3: return 0.0;
619  }
620  return 0.0;
621 }
Vector bb_max
Definition: ex9.cpp:67
void SetPrecision(int prec)
Set the precision (number of digits) used for the text output of doubles.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:563
Conjugate gradient method.
Definition: solvers.hpp:465
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Data type for scaled Jacobi-type smoother of sparse matrix.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void SetDataFormat(VTKFormat fmt)
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:711
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
Helper class for ParaView visualization data.
Base abstract class for first order time dependent operators.
Definition: operator.hpp:285
AssemblyLevel GetAssemblyLevel() const
Returns the assembly level.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition: mesh.cpp:129
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:102
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]...
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Abstract class for solving systems of ODEs: dx/dt = f(x,t)
Definition: ode.hpp:22
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
virtual void SetTime(const double t_)
Set the current time.
Definition: operator.hpp:329
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:655
FDualNumber< tbase > exp(const FDualNumber< tbase > &f)
exp([dual number])
Definition: fdual.hpp:499
Backward Euler ODE solver. L-stable.
Definition: ode.hpp:391
Data collection with Sidre routines following the Conduit mesh blueprint specification.
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3619
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:46
constexpr char vishost[]
virtual void Save()
Save the collection to disk.
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
constexpr int visport
Data collection with VisIt I/O routines.
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:5312
void SetHighOrderOutput(bool high_order_output_)
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:83
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetTime(double t)
Set physical time (for time-dependent simulations)
Vector bb_min
Definition: ex9.cpp:67
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
A general vector function coefficient.
The classical explicit forth-order Runge-Kutta method, RK4.
Definition: ode.hpp:162
void SetAbsTol(double atol)
Definition: solvers.hpp:199
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Definition: ex9.cpp:484
void SetRelTol(double rtol)
Definition: solvers.hpp:198
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
void velocity_function(const Vector &x, Vector &v)
Definition: ex9.cpp:501
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
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)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
int problem
Definition: ex15.cpp:62
Third-order, strong stability preserving (SSP) Runge-Kutta method.
Definition: ode.hpp:149
GMRES method.
Definition: solvers.hpp:497
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:128
Implicit midpoint method. A-stable, not L-stable.
Definition: ode.hpp:404
FE_Evolution(FiniteElementSpace &vfes_, Operator &A_, SparseMatrix &Aflux_)
Definition: ex18.hpp:102
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:272
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1646
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
int dim
Definition: ex24.cpp:53
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2395
virtual ~FE_Evolution()
Definition: ex18.hpp:43
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
double u0_function(const Vector &x)
Definition: ex9.cpp:558
RefCoord t[3]
void SetLevelsOfDetail(int levels_of_detail_)
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
virtual void Save() override
RefCoord s[3]
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix.
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Base class for solvers.
Definition: operator.hpp:651
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
The classical forward Euler method.
Definition: ode.hpp:116
Abstract operator.
Definition: operator.hpp:24
virtual void Init(TimeDependentOperator &f_)
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:18
int main()
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
double inflow_function(const Vector &x)
Definition: ex9.cpp:611
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:284
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
alpha (q . grad u, v)