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