MFEM  v4.6.0
Finite element discretization library
ex22.cpp
Go to the documentation of this file.
1 // MFEM Example 22
2 //
3 // Compile with: make ex22
4 //
5 // Sample runs: ex22 -m ../data/inline-segment.mesh -o 3
6 // ex22 -m ../data/inline-tri.mesh -o 3
7 // ex22 -m ../data/inline-quad.mesh -o 3
8 // ex22 -m ../data/inline-quad.mesh -o 3 -p 1
9 // ex22 -m ../data/inline-quad.mesh -o 3 -p 1 -pa
10 // ex22 -m ../data/inline-quad.mesh -o 3 -p 2
11 // ex22 -m ../data/inline-tet.mesh -o 2
12 // ex22 -m ../data/inline-hex.mesh -o 2
13 // ex22 -m ../data/inline-hex.mesh -o 2 -p 1
14 // ex22 -m ../data/inline-hex.mesh -o 2 -p 2
15 // ex22 -m ../data/inline-hex.mesh -o 2 -p 2 -pa
16 // ex22 -m ../data/inline-wedge.mesh -o 1
17 // ex22 -m ../data/inline-pyramid.mesh -o 1
18 // ex22 -m ../data/star.mesh -r 1 -o 2 -sigma 10.0
19 //
20 // Device sample runs:
21 // ex22 -m ../data/inline-quad.mesh -o 3 -p 1 -pa -d cuda
22 // ex22 -m ../data/inline-hex.mesh -o 2 -p 2 -pa -d cuda
23 // ex22 -m ../data/star.mesh -r 1 -o 2 -sigma 10.0 -pa -d cuda
24 //
25 // Description: This example code demonstrates the use of MFEM to define and
26 // solve simple complex-valued linear systems. It implements three
27 // variants of a damped harmonic oscillator:
28 //
29 // 1) A scalar H1 field
30 // -Div(a Grad u) - omega^2 b u + i omega c u = 0
31 //
32 // 2) A vector H(Curl) field
33 // Curl(a Curl u) - omega^2 b u + i omega c u = 0
34 //
35 // 3) A vector H(Div) field
36 // -Grad(a Div u) - omega^2 b u + i omega c u = 0
37 //
38 // In each case the field is driven by a forced oscillation, with
39 // angular frequency omega, imposed at the boundary or a portion
40 // of the boundary.
41 //
42 // In electromagnetics, the coefficients are typically named the
43 // permeability, mu = 1/a, permittivity, epsilon = b, and
44 // conductivity, sigma = c. The user can specify these constants
45 // using either set of names.
46 //
47 // The example also demonstrates how to display a time-varying
48 // solution as a sequence of fields sent to a single GLVis socket.
49 //
50 // We recommend viewing examples 1, 3 and 4 before viewing this
51 // example.
52 
53 #include "mfem.hpp"
54 #include <fstream>
55 #include <iostream>
56 
57 using namespace std;
58 using namespace mfem;
59 
60 static double mu_ = 1.0;
61 static double epsilon_ = 1.0;
62 static double sigma_ = 20.0;
63 static double omega_ = 10.0;
64 
65 double u0_real_exact(const Vector &);
66 double u0_imag_exact(const Vector &);
67 
68 void u1_real_exact(const Vector &, Vector &);
69 void u1_imag_exact(const Vector &, Vector &);
70 
71 void u2_real_exact(const Vector &, Vector &);
72 void u2_imag_exact(const Vector &, Vector &);
73 
74 bool check_for_inline_mesh(const char * mesh_file);
75 
76 int main(int argc, char *argv[])
77 {
78  // 1. Parse command-line options.
79  const char *mesh_file = "../data/inline-quad.mesh";
80  int ref_levels = 0;
81  int order = 1;
82  int prob = 0;
83  double freq = -1.0;
84  double a_coef = 0.0;
85  bool visualization = 1;
86  bool herm_conv = true;
87  bool exact_sol = true;
88  bool pa = false;
89  const char *device_config = "cpu";
90 
91  OptionsParser args(argc, argv);
92  args.AddOption(&mesh_file, "-m", "--mesh",
93  "Mesh file to use.");
94  args.AddOption(&ref_levels, "-r", "--refine",
95  "Number of times to refine the mesh uniformly.");
96  args.AddOption(&order, "-o", "--order",
97  "Finite element order (polynomial degree).");
98  args.AddOption(&prob, "-p", "--problem-type",
99  "Choose between 0: H_1, 1: H(Curl), or 2: H(Div) "
100  "damped harmonic oscillator.");
101  args.AddOption(&a_coef, "-a", "--stiffness-coef",
102  "Stiffness coefficient (spring constant or 1/mu).");
103  args.AddOption(&epsilon_, "-b", "--mass-coef",
104  "Mass coefficient (or epsilon).");
105  args.AddOption(&sigma_, "-c", "--damping-coef",
106  "Damping coefficient (or sigma).");
107  args.AddOption(&mu_, "-mu", "--permeability",
108  "Permeability of free space (or 1/(spring constant)).");
109  args.AddOption(&epsilon_, "-eps", "--permittivity",
110  "Permittivity of free space (or mass constant).");
111  args.AddOption(&sigma_, "-sigma", "--conductivity",
112  "Conductivity (or damping constant).");
113  args.AddOption(&freq, "-f", "--frequency",
114  "Frequency (in Hz).");
115  args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
116  "--no-hermitian", "Use convention for Hermitian operators.");
117  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
118  "--no-visualization",
119  "Enable or disable GLVis visualization.");
120  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
121  "--no-partial-assembly", "Enable Partial Assembly.");
122  args.AddOption(&device_config, "-d", "--device",
123  "Device configuration string, see Device::Configure().");
124  args.Parse();
125  if (!args.Good())
126  {
127  args.PrintUsage(cout);
128  return 1;
129  }
130  args.PrintOptions(cout);
131 
132  MFEM_VERIFY(prob >= 0 && prob <=2,
133  "Unrecognized problem type: " << prob);
134 
135  if ( a_coef != 0.0 )
136  {
137  mu_ = 1.0 / a_coef;
138  }
139  if ( freq > 0.0 )
140  {
141  omega_ = 2.0 * M_PI * freq;
142  }
143 
144  exact_sol = check_for_inline_mesh(mesh_file);
145  if (exact_sol)
146  {
147  cout << "Identified a mesh with known exact solution" << endl;
148  }
149 
151  herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
152 
153  // 2. Enable hardware devices such as GPUs, and programming models such as
154  // CUDA, OCCA, RAJA and OpenMP based on command line options.
155  Device device(device_config);
156  device.Print();
157 
158  // 3. Read the mesh from the given mesh file. We can handle triangular,
159  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes
160  // with the same code.
161  Mesh *mesh = new Mesh(mesh_file, 1, 1);
162  int dim = mesh->Dimension();
163 
164  // 4. Refine the mesh to increase resolution. In this example we do
165  // 'ref_levels' of uniform refinement where the user specifies
166  // the number of levels with the '-r' option.
167  for (int l = 0; l < ref_levels; l++)
168  {
169  mesh->UniformRefinement();
170  }
171 
172  // 5. Define a finite element space on the mesh. Here we use continuous
173  // Lagrange, Nedelec, or Raviart-Thomas finite elements of the specified
174  // order.
175  if (dim == 1 && prob != 0 )
176  {
177  cout << "Switching to problem type 0, H1 basis functions, "
178  << "for 1 dimensional mesh." << endl;
179  prob = 0;
180  }
181 
182  FiniteElementCollection *fec = NULL;
183  switch (prob)
184  {
185  case 0: fec = new H1_FECollection(order, dim); break;
186  case 1: fec = new ND_FECollection(order, dim); break;
187  case 2: fec = new RT_FECollection(order - 1, dim); break;
188  default: break; // This should be unreachable
189  }
190  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
191  cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
192  << endl;
193 
194  // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
195  // In this example, the boundary conditions are defined based on the type
196  // of mesh and the problem type.
197  Array<int> ess_tdof_list;
198  Array<int> ess_bdr;
199  if (mesh->bdr_attributes.Size())
200  {
201  ess_bdr.SetSize(mesh->bdr_attributes.Max());
202  ess_bdr = 1;
203  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
204  }
205 
206  // 7. Set up the linear form b(.) which corresponds to the right-hand side of
207  // the FEM linear system.
208  ComplexLinearForm b(fespace, conv);
209  b.Vector::operator=(0.0);
210 
211  // 8. Define the solution vector u as a complex finite element grid function
212  // corresponding to fespace. Initialize u with initial guess of 1+0i or
213  // the exact solution if it is known.
214  ComplexGridFunction u(fespace);
215  ComplexGridFunction * u_exact = NULL;
216  if (exact_sol) { u_exact = new ComplexGridFunction(fespace); }
217 
224 
225  ConstantCoefficient zeroCoef(0.0);
226  ConstantCoefficient oneCoef(1.0);
227 
228  Vector zeroVec(dim); zeroVec = 0.0;
229  Vector oneVec(dim); oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
230  VectorConstantCoefficient zeroVecCoef(zeroVec);
231  VectorConstantCoefficient oneVecCoef(oneVec);
232 
233  switch (prob)
234  {
235  case 0:
236  if (exact_sol)
237  {
238  u.ProjectBdrCoefficient(u0_r, u0_i, ess_bdr);
239  u_exact->ProjectCoefficient(u0_r, u0_i);
240  }
241  else
242  {
243  u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
244  }
245  break;
246  case 1:
247  if (exact_sol)
248  {
249  u.ProjectBdrCoefficientTangent(u1_r, u1_i, ess_bdr);
250  u_exact->ProjectCoefficient(u1_r, u1_i);
251  }
252  else
253  {
254  u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
255  }
256  break;
257  case 2:
258  if (exact_sol)
259  {
260  u.ProjectBdrCoefficientNormal(u2_r, u2_i, ess_bdr);
261  u_exact->ProjectCoefficient(u2_r, u2_i);
262  }
263  else
264  {
265  u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
266  }
267  break;
268  default: break; // This should be unreachable
269  }
270 
271  if (visualization && exact_sol)
272  {
273  char vishost[] = "localhost";
274  int visport = 19916;
275  socketstream sol_sock_r(vishost, visport);
276  socketstream sol_sock_i(vishost, visport);
277  sol_sock_r.precision(8);
278  sol_sock_i.precision(8);
279  sol_sock_r << "solution\n" << *mesh << u_exact->real()
280  << "window_title 'Exact: Real Part'" << flush;
281  sol_sock_i << "solution\n" << *mesh << u_exact->imag()
282  << "window_title 'Exact: Imaginary Part'" << flush;
283  }
284 
285  // 9. Set up the sesquilinear form a(.,.) on the finite element space
286  // corresponding to the damped harmonic oscillator operator of the
287  // appropriate type:
288  //
289  // 0) A scalar H1 field
290  // -Div(a Grad) - omega^2 b + i omega c
291  //
292  // 1) A vector H(Curl) field
293  // Curl(a Curl) - omega^2 b + i omega c
294  //
295  // 2) A vector H(Div) field
296  // -Grad(a Div) - omega^2 b + i omega c
297  //
298  ConstantCoefficient stiffnessCoef(1.0/mu_);
299  ConstantCoefficient massCoef(-omega_ * omega_ * epsilon_);
300  ConstantCoefficient lossCoef(omega_ * sigma_);
301  ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
302 
303  SesquilinearForm *a = new SesquilinearForm(fespace, conv);
304  if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
305  switch (prob)
306  {
307  case 0:
308  a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
309  NULL);
310  a->AddDomainIntegrator(new MassIntegrator(massCoef),
311  new MassIntegrator(lossCoef));
312  break;
313  case 1:
314  a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
315  NULL);
316  a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
317  new VectorFEMassIntegrator(lossCoef));
318  break;
319  case 2:
320  a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
321  NULL);
322  a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
323  new VectorFEMassIntegrator(lossCoef));
324  break;
325  default: break; // This should be unreachable
326  }
327 
328  // 9a. Set up the bilinear form for the preconditioner corresponding to the
329  // appropriate operator
330  //
331  // 0) A scalar H1 field
332  // -Div(a Grad) - omega^2 b + omega c
333  //
334  // 1) A vector H(Curl) field
335  // Curl(a Curl) + omega^2 b + omega c
336  //
337  // 2) A vector H(Div) field
338  // -Grad(a Div) - omega^2 b + omega c
339  //
340  BilinearForm *pcOp = new BilinearForm(fespace);
341  if (pa) { pcOp->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
342 
343  switch (prob)
344  {
345  case 0:
346  pcOp->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef));
347  pcOp->AddDomainIntegrator(new MassIntegrator(massCoef));
348  pcOp->AddDomainIntegrator(new MassIntegrator(lossCoef));
349  break;
350  case 1:
351  pcOp->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef));
352  pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(negMassCoef));
353  pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
354  break;
355  case 2:
356  pcOp->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef));
357  pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef));
358  pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
359  break;
360  default: break; // This should be unreachable
361  }
362 
363  // 10. Assemble the form and the corresponding linear system, applying any
364  // necessary transformations such as: assembly, eliminating boundary
365  // conditions, conforming constraints for non-conforming AMR, etc.
366  a->Assemble();
367  pcOp->Assemble();
368 
369  OperatorHandle A;
370  Vector B, U;
371 
372  a->FormLinearSystem(ess_tdof_list, u, b, A, U, B);
373 
374  cout << "Size of linear system: " << A->Width() << endl << endl;
375 
376  // 11. Define and apply a GMRES solver for AU=B with a block diagonal
377  // preconditioner based on the appropriate sparse smoother.
378  {
379  Array<int> blockOffsets;
380  blockOffsets.SetSize(3);
381  blockOffsets[0] = 0;
382  blockOffsets[1] = A->Height() / 2;
383  blockOffsets[2] = A->Height() / 2;
384  blockOffsets.PartialSum();
385 
386  BlockDiagonalPreconditioner BDP(blockOffsets);
387 
388  Operator * pc_r = NULL;
389  Operator * pc_i = NULL;
390 
391  if (pa)
392  {
393  pc_r = new OperatorJacobiSmoother(*pcOp, ess_tdof_list);
394  }
395  else
396  {
397  OperatorHandle PCOp;
399  pcOp->FormSystemMatrix(ess_tdof_list, PCOp);
400  switch (prob)
401  {
402  case 0:
403  pc_r = new DSmoother(*PCOp.As<SparseMatrix>());
404  break;
405  case 1:
406  pc_r = new GSSmoother(*PCOp.As<SparseMatrix>());
407  break;
408  case 2:
409  pc_r = new DSmoother(*PCOp.As<SparseMatrix>());
410  break;
411  default:
412  break; // This should be unreachable
413  }
414  }
415  double s = (prob != 1) ? 1.0 : -1.0;
416  pc_i = new ScaledOperator(pc_r,
417  (conv == ComplexOperator::HERMITIAN) ?
418  s:-s);
419 
420  BDP.SetDiagonalBlock(0, pc_r);
421  BDP.SetDiagonalBlock(1, pc_i);
422  BDP.owns_blocks = 1;
423 
424  GMRESSolver gmres;
425  gmres.SetPreconditioner(BDP);
426  gmres.SetOperator(*A.Ptr());
427  gmres.SetRelTol(1e-12);
428  gmres.SetMaxIter(1000);
429  gmres.SetPrintLevel(1);
430  gmres.Mult(B, U);
431  }
432 
433  // 12. Recover the solution as a finite element grid function and compute the
434  // errors if the exact solution is known.
435  a->RecoverFEMSolution(U, b, u);
436 
437  if (exact_sol)
438  {
439  double err_r = -1.0;
440  double err_i = -1.0;
441 
442  switch (prob)
443  {
444  case 0:
445  err_r = u.real().ComputeL2Error(u0_r);
446  err_i = u.imag().ComputeL2Error(u0_i);
447  break;
448  case 1:
449  err_r = u.real().ComputeL2Error(u1_r);
450  err_i = u.imag().ComputeL2Error(u1_i);
451  break;
452  case 2:
453  err_r = u.real().ComputeL2Error(u2_r);
454  err_i = u.imag().ComputeL2Error(u2_i);
455  break;
456  default: break; // This should be unreachable
457  }
458 
459  cout << endl;
460  cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
461  cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
462  cout << endl;
463  }
464 
465  // 13. Save the refined mesh and the solution. This output can be viewed
466  // later using GLVis: "glvis -m mesh -g sol".
467  {
468  ofstream mesh_ofs("refined.mesh");
469  mesh_ofs.precision(8);
470  mesh->Print(mesh_ofs);
471 
472  ofstream sol_r_ofs("sol_r.gf");
473  ofstream sol_i_ofs("sol_i.gf");
474  sol_r_ofs.precision(8);
475  sol_i_ofs.precision(8);
476  u.real().Save(sol_r_ofs);
477  u.imag().Save(sol_i_ofs);
478  }
479 
480  // 14. Send the solution by socket to a GLVis server.
481  if (visualization)
482  {
483  char vishost[] = "localhost";
484  int visport = 19916;
485  socketstream sol_sock_r(vishost, visport);
486  socketstream sol_sock_i(vishost, visport);
487  sol_sock_r.precision(8);
488  sol_sock_i.precision(8);
489  sol_sock_r << "solution\n" << *mesh << u.real()
490  << "window_title 'Solution: Real Part'" << flush;
491  sol_sock_i << "solution\n" << *mesh << u.imag()
492  << "window_title 'Solution: Imaginary Part'" << flush;
493  }
494  if (visualization && exact_sol)
495  {
496  *u_exact -= u;
497 
498  char vishost[] = "localhost";
499  int visport = 19916;
500  socketstream sol_sock_r(vishost, visport);
501  socketstream sol_sock_i(vishost, visport);
502  sol_sock_r.precision(8);
503  sol_sock_i.precision(8);
504  sol_sock_r << "solution\n" << *mesh << u_exact->real()
505  << "window_title 'Error: Real Part'" << flush;
506  sol_sock_i << "solution\n" << *mesh << u_exact->imag()
507  << "window_title 'Error: Imaginary Part'" << flush;
508  }
509  if (visualization)
510  {
511  GridFunction u_t(fespace);
512  u_t = u.real();
513  char vishost[] = "localhost";
514  int visport = 19916;
515  socketstream sol_sock(vishost, visport);
516  sol_sock.precision(8);
517  sol_sock << "solution\n" << *mesh << u_t
518  << "window_title 'Harmonic Solution (t = 0.0 T)'"
519  << "pause\n" << flush;
520 
521  cout << "GLVis visualization paused."
522  << " Press space (in the GLVis window) to resume it.\n";
523  int num_frames = 32;
524  int i = 0;
525  while (sol_sock)
526  {
527  double t = (double)(i % num_frames) / num_frames;
528  ostringstream oss;
529  oss << "Harmonic Solution (t = " << t << " T)";
530 
531  add(cos( 2.0 * M_PI * t), u.real(),
532  sin(-2.0 * M_PI * t), u.imag(), u_t);
533  sol_sock << "solution\n" << *mesh << u_t
534  << "window_title '" << oss.str() << "'" << flush;
535  i++;
536  }
537  }
538 
539  // 15. Free the used memory.
540  delete a;
541  delete u_exact;
542  delete pcOp;
543  delete fespace;
544  delete fec;
545  delete mesh;
546 
547  return 0;
548 }
549 
550 bool check_for_inline_mesh(const char * mesh_file)
551 {
552  string file(mesh_file);
553  size_t p0 = file.find_last_of("/");
554  string s0 = file.substr((p0==string::npos)?0:(p0+1),7);
555  return s0 == "inline-";
556 }
557 
558 complex<double> u0_exact(const Vector &x)
559 {
560  int dim = x.Size();
561  complex<double> i(0.0, 1.0);
562  complex<double> alpha = (epsilon_ * omega_ - i * sigma_);
563  complex<double> kappa = std::sqrt(mu_ * omega_* alpha);
564  return std::exp(-i * kappa * x[dim - 1]);
565 }
566 
567 double u0_real_exact(const Vector &x)
568 {
569  return u0_exact(x).real();
570 }
571 
572 double u0_imag_exact(const Vector &x)
573 {
574  return u0_exact(x).imag();
575 }
576 
577 void u1_real_exact(const Vector &x, Vector &v)
578 {
579  int dim = x.Size();
580  v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
581 }
582 
583 void u1_imag_exact(const Vector &x, Vector &v)
584 {
585  int dim = x.Size();
586  v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
587 }
588 
589 void u2_real_exact(const Vector &x, Vector &v)
590 {
591  int dim = x.Size();
592  v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
593 }
594 
595 void u2_imag_exact(const Vector &x, Vector &v)
596 {
597  int dim = x.Size();
598  v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
599 }
int visport
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
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.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void u2_imag_exact(const Vector &, Vector &)
Definition: ex22.cpp:595
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Integrator for (curl u, curl v) for Nedelec elements.
Vector coefficient that is constant in space and time.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:980
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:587
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.
double kappa
Definition: ex24.cpp:54
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
complex< double > u0_exact(const Vector &x)
Definition: ex22.cpp:558
STL namespace.
(Q div u, div v) for RT elements
void u1_imag_exact(const Vector &, Vector &)
Definition: ex22.cpp:583
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:317
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
Data type for Gauss-Seidel smoother of sparse matrix.
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
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
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
Data type sparse matrix.
Definition: sparsemat.hpp:50
char vishost[]
int main(int argc, char *argv[])
Definition: ex22.cpp:76
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
double u0_imag_exact(const Vector &)
Definition: ex22.cpp:572
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Set the diagonal value to one.
Definition: operator.hpp:50
prob_type prob
Definition: ex25.cpp:154
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
A general vector function coefficient.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void u_exact(const Vector &x, Vector &u)
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Scaled Operator B: x -> a A(x).
Definition: operator.hpp:726
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
GMRES method.
Definition: solvers.hpp:525
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
double a
Definition: lissajous.cpp:41
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
int dim
Definition: ex24.cpp:53
double freq
Definition: ex24.cpp:54
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:454
Vector data type.
Definition: vector.hpp:58
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:2011
void u1_real_exact(const Vector &, Vector &)
Definition: ex22.cpp:577
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
RefCoord s[3]
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Abstract operator.
Definition: operator.hpp:24
bool check_for_inline_mesh(const char *mesh_file)
Definition: ex22.cpp:550
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
void u2_real_exact(const Vector &, Vector &)
Definition: ex22.cpp:589
double u0_real_exact(const Vector &)
Definition: ex22.cpp:567