MFEM  v4.6.0
Finite element discretization library
ex22p.cpp
Go to the documentation of this file.
1 // MFEM Example 22 - Parallel Version
2 //
3 // Compile with: make ex22p
4 //
5 // Sample runs: mpirun -np 4 ex22p -m ../data/inline-segment.mesh -o 3
6 // mpirun -np 4 ex22p -m ../data/inline-tri.mesh -o 3
7 // mpirun -np 4 ex22p -m ../data/inline-quad.mesh -o 3
8 // mpirun -np 4 ex22p -m ../data/inline-quad.mesh -o 3 -p 1
9 // mpirun -np 4 ex22p -m ../data/inline-quad.mesh -o 3 -p 2
10 // mpirun -np 4 ex22p -m ../data/inline-quad.mesh -o 1 -p 1 -pa
11 // mpirun -np 4 ex22p -m ../data/inline-tet.mesh -o 2
12 // mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 2
13 // mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 2 -p 1
14 // mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 2 -p 2
15 // mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 1 -p 2 -pa
16 // mpirun -np 4 ex22p -m ../data/inline-wedge.mesh -o 1
17 // mpirun -np 4 ex22p -m ../data/inline-pyramid.mesh -o 1
18 // mpirun -np 4 ex22p -m ../data/star.mesh -o 2 -sigma 10.0
19 //
20 // Device sample runs:
21 // mpirun -np 4 ex22p -m ../data/inline-quad.mesh -o 1 -p 1 -pa -d cuda
22 // mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 1 -p 2 -pa -d cuda
23 // mpirun -np 4 ex22p -m ../data/star.mesh -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. Initialize MPI and HYPRE.
79  Mpi::Init(argc, argv);
80  int num_procs = Mpi::WorldSize();
81  int myid = Mpi::WorldRank();
82  Hypre::Init();
83 
84  // 2. Parse command-line options.
85  const char *mesh_file = "../data/inline-quad.mesh";
86  int ser_ref_levels = 1;
87  int par_ref_levels = 1;
88  int order = 1;
89  int prob = 0;
90  double freq = -1.0;
91  double a_coef = 0.0;
92  bool visualization = 1;
93  bool herm_conv = true;
94  bool exact_sol = true;
95  bool pa = false;
96  const char *device_config = "cpu";
97 
98  OptionsParser args(argc, argv);
99  args.AddOption(&mesh_file, "-m", "--mesh",
100  "Mesh file to use.");
101  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
102  "Number of times to refine the mesh uniformly in serial.");
103  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
104  "Number of times to refine the mesh uniformly in parallel.");
105  args.AddOption(&order, "-o", "--order",
106  "Finite element order (polynomial degree).");
107  args.AddOption(&prob, "-p", "--problem-type",
108  "Choose between 0: H_1, 1: H(Curl), or 2: H(Div) "
109  "damped harmonic oscillator.");
110  args.AddOption(&a_coef, "-a", "--stiffness-coef",
111  "Stiffness coefficient (spring constant or 1/mu).");
112  args.AddOption(&epsilon_, "-b", "--mass-coef",
113  "Mass coefficient (or epsilon).");
114  args.AddOption(&sigma_, "-c", "--damping-coef",
115  "Damping coefficient (or sigma).");
116  args.AddOption(&mu_, "-mu", "--permeability",
117  "Permeability of free space (or 1/(spring constant)).");
118  args.AddOption(&epsilon_, "-eps", "--permittivity",
119  "Permittivity of free space (or mass constant).");
120  args.AddOption(&sigma_, "-sigma", "--conductivity",
121  "Conductivity (or damping constant).");
122  args.AddOption(&freq, "-f", "--frequency",
123  "Frequency (in Hz).");
124  args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
125  "--no-hermitian", "Use convention for Hermitian operators.");
126  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
127  "--no-visualization",
128  "Enable or disable GLVis visualization.");
129  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
130  "--no-partial-assembly", "Enable Partial Assembly.");
131  args.AddOption(&device_config, "-d", "--device",
132  "Device configuration string, see Device::Configure().");
133  args.Parse();
134  if (!args.Good())
135  {
136  if (myid == 0)
137  {
138  args.PrintUsage(cout);
139  }
140  return 1;
141  }
142  if (myid == 0)
143  {
144  args.PrintOptions(cout);
145  }
146 
147  MFEM_VERIFY(prob >= 0 && prob <=2,
148  "Unrecognized problem type: " << prob);
149 
150  if ( a_coef != 0.0 )
151  {
152  mu_ = 1.0 / a_coef;
153  }
154  if ( freq > 0.0 )
155  {
156  omega_ = 2.0 * M_PI * freq;
157  }
158 
159  exact_sol = check_for_inline_mesh(mesh_file);
160  if (myid == 0 && exact_sol)
161  {
162  cout << "Identified a mesh with known exact solution" << endl;
163  }
164 
166  herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
167 
168  // 3. Enable hardware devices such as GPUs, and programming models such as
169  // CUDA, OCCA, RAJA and OpenMP based on command line options.
170  Device device(device_config);
171  if (myid == 0) { device.Print(); }
172 
173  // 4. Read the (serial) mesh from the given mesh file on all processors. We
174  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
175  // and volume meshes with the same code.
176  Mesh *mesh = new Mesh(mesh_file, 1, 1);
177  int dim = mesh->Dimension();
178 
179  // 5. Refine the serial mesh on all processors to increase the resolution.
180  for (int l = 0; l < ser_ref_levels; l++)
181  {
182  mesh->UniformRefinement();
183  }
184 
185  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
186  // this mesh further in parallel to increase the resolution. Once the
187  // parallel mesh is defined, the serial mesh can be deleted.
188  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
189  delete mesh;
190  for (int l = 0; l < par_ref_levels; l++)
191  {
192  pmesh->UniformRefinement();
193  }
194 
195  // 7. Define a parallel finite element space on the parallel mesh. Here we
196  // use continuous Lagrange, Nedelec, or Raviart-Thomas finite elements of
197  // the specified order.
198  if (dim == 1 && prob != 0 )
199  {
200  if (myid == 0)
201  {
202  cout << "Switching to problem type 0, H1 basis functions, "
203  << "for 1 dimensional mesh." << endl;
204  }
205  prob = 0;
206  }
207 
208  FiniteElementCollection *fec = NULL;
209  switch (prob)
210  {
211  case 0: fec = new H1_FECollection(order, dim); break;
212  case 1: fec = new ND_FECollection(order, dim); break;
213  case 2: fec = new RT_FECollection(order - 1, dim); break;
214  default: break; // This should be unreachable
215  }
216  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
217  HYPRE_BigInt size = fespace->GlobalTrueVSize();
218  if (myid == 0)
219  {
220  cout << "Number of finite element unknowns: " << size << endl;
221  }
222 
223  // 8. Determine the list of true (i.e. parallel conforming) essential
224  // boundary dofs. In this example, the boundary conditions are defined
225  // based on the type of mesh and the problem type.
226  Array<int> ess_tdof_list;
227  Array<int> ess_bdr;
228  if (pmesh->bdr_attributes.Size())
229  {
230  ess_bdr.SetSize(pmesh->bdr_attributes.Max());
231  ess_bdr = 1;
232  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
233  }
234 
235  // 9. Set up the parallel linear form b(.) which corresponds to the
236  // right-hand side of the FEM linear system.
237  ParComplexLinearForm b(fespace, conv);
238  b.Vector::operator=(0.0);
239 
240  // 10. Define the solution vector u as a parallel complex finite element grid
241  // function corresponding to fespace. Initialize u with initial guess of
242  // 1+0i or the exact solution if it is known.
243  ParComplexGridFunction u(fespace);
245  if (exact_sol) { u_exact = new ParComplexGridFunction(fespace); }
246 
253 
254  ConstantCoefficient zeroCoef(0.0);
255  ConstantCoefficient oneCoef(1.0);
256 
257  Vector zeroVec(dim); zeroVec = 0.0;
258  Vector oneVec(dim); oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
259  VectorConstantCoefficient zeroVecCoef(zeroVec);
260  VectorConstantCoefficient oneVecCoef(oneVec);
261 
262  switch (prob)
263  {
264  case 0:
265  if (exact_sol)
266  {
267  u.ProjectBdrCoefficient(u0_r, u0_i, ess_bdr);
268  u_exact->ProjectCoefficient(u0_r, u0_i);
269  }
270  else
271  {
272  u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
273  }
274  break;
275  case 1:
276  if (exact_sol)
277  {
278  u.ProjectBdrCoefficientTangent(u1_r, u1_i, ess_bdr);
279  u_exact->ProjectCoefficient(u1_r, u1_i);
280  }
281  else
282  {
283  u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
284  }
285  break;
286  case 2:
287  if (exact_sol)
288  {
289  u.ProjectBdrCoefficientNormal(u2_r, u2_i, ess_bdr);
290  u_exact->ProjectCoefficient(u2_r, u2_i);
291  }
292  else
293  {
294  u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
295  }
296  break;
297  default: break; // This should be unreachable
298  }
299 
300  if (visualization && exact_sol)
301  {
302  char vishost[] = "localhost";
303  int visport = 19916;
304  socketstream sol_sock_r(vishost, visport);
305  socketstream sol_sock_i(vishost, visport);
306  sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
307  sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
308  sol_sock_r.precision(8);
309  sol_sock_i.precision(8);
310  sol_sock_r << "solution\n" << *pmesh << u_exact->real()
311  << "window_title 'Exact: Real Part'" << flush;
312  sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
313  << "window_title 'Exact: Imaginary Part'" << flush;
314  }
315 
316  // 11. Set up the parallel sesquilinear form a(.,.) on the finite element
317  // space corresponding to the damped harmonic oscillator operator of the
318  // appropriate type:
319  //
320  // 0) A scalar H1 field
321  // -Div(a Grad) - omega^2 b + i omega c
322  //
323  // 1) A vector H(Curl) field
324  // Curl(a Curl) - omega^2 b + i omega c
325  //
326  // 2) A vector H(Div) field
327  // -Grad(a Div) - omega^2 b + i omega c
328  //
329  ConstantCoefficient stiffnessCoef(1.0/mu_);
330  ConstantCoefficient massCoef(-omega_ * omega_ * epsilon_);
331  ConstantCoefficient lossCoef(omega_ * sigma_);
332  ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
333 
334  ParSesquilinearForm *a = new ParSesquilinearForm(fespace, conv);
335  if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
336  switch (prob)
337  {
338  case 0:
339  a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
340  NULL);
341  a->AddDomainIntegrator(new MassIntegrator(massCoef),
342  new MassIntegrator(lossCoef));
343  break;
344  case 1:
345  a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
346  NULL);
347  a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
348  new VectorFEMassIntegrator(lossCoef));
349  break;
350  case 2:
351  a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
352  NULL);
353  a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
354  new VectorFEMassIntegrator(lossCoef));
355  break;
356  default: break; // This should be unreachable
357  }
358 
359  // 11a. Set up the parallel bilinear form for the preconditioner
360  // corresponding to the appropriate operator
361  //
362  // 0) A scalar H1 field
363  // -Div(a Grad) - omega^2 b + omega c
364  //
365  // 1) A vector H(Curl) field
366  // Curl(a Curl) + omega^2 b + omega c
367  //
368  // 2) A vector H(Div) field
369  // -Grad(a Div) - omega^2 b + omega c
370  //
371  ParBilinearForm *pcOp = new ParBilinearForm(fespace);
372  if (pa) { pcOp->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
373  switch (prob)
374  {
375  case 0:
376  pcOp->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef));
377  pcOp->AddDomainIntegrator(new MassIntegrator(massCoef));
378  pcOp->AddDomainIntegrator(new MassIntegrator(lossCoef));
379  break;
380  case 1:
381  pcOp->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef));
382  pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(negMassCoef));
383  pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
384  break;
385  case 2:
386  pcOp->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef));
387  pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef));
388  pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
389  break;
390  default: break; // This should be unreachable
391  }
392 
393  // 12. Assemble the parallel bilinear form and the corresponding linear
394  // system, applying any necessary transformations such as: parallel
395  // assembly, eliminating boundary conditions, applying conforming
396  // constraints for non-conforming AMR, etc.
397  a->Assemble();
398  pcOp->Assemble();
399 
400  OperatorHandle A;
401  Vector B, U;
402 
403  a->FormLinearSystem(ess_tdof_list, u, b, A, U, B);
404 
405  if (myid == 0)
406  {
407  cout << "Size of linear system: "
408  << 2 * fespace->GlobalTrueVSize() << endl << endl;
409  }
410 
411  // 13. Define and apply a parallel FGMRES solver for AU=B with a block
412  // diagonal preconditioner based on the appropriate multigrid
413  // preconditioner from hypre.
414  {
415  Array<int> blockTrueOffsets;
416  blockTrueOffsets.SetSize(3);
417  blockTrueOffsets[0] = 0;
418  blockTrueOffsets[1] = A->Height() / 2;
419  blockTrueOffsets[2] = A->Height() / 2;
420  blockTrueOffsets.PartialSum();
421 
422  BlockDiagonalPreconditioner BDP(blockTrueOffsets);
423 
424  Operator * pc_r = NULL;
425  Operator * pc_i = NULL;
426 
427  if (pa)
428  {
429  pc_r = new OperatorJacobiSmoother(*pcOp, ess_tdof_list);
430  }
431  else
432  {
433  OperatorHandle PCOp;
434  pcOp->FormSystemMatrix(ess_tdof_list, PCOp);
435  switch (prob)
436  {
437  case 0:
438  pc_r = new HypreBoomerAMG(*PCOp.As<HypreParMatrix>());
439  break;
440  case 1:
441  pc_r = new HypreAMS(*PCOp.As<HypreParMatrix>(), fespace);
442  break;
443  case 2:
444  if (dim == 2 )
445  {
446  pc_r = new HypreAMS(*PCOp.As<HypreParMatrix>(), fespace);
447  }
448  else
449  {
450  pc_r = new HypreADS(*PCOp.As<HypreParMatrix>(), fespace);
451  }
452  break;
453  default: break; // This should be unreachable
454  }
455  }
456  pc_i = new ScaledOperator(pc_r,
457  (conv == ComplexOperator::HERMITIAN) ?
458  -1.0:1.0);
459 
460  BDP.SetDiagonalBlock(0, pc_r);
461  BDP.SetDiagonalBlock(1, pc_i);
462  BDP.owns_blocks = 1;
463 
464  FGMRESSolver fgmres(MPI_COMM_WORLD);
465  fgmres.SetPreconditioner(BDP);
466  fgmres.SetOperator(*A.Ptr());
467  fgmres.SetRelTol(1e-12);
468  fgmres.SetMaxIter(1000);
469  fgmres.SetPrintLevel(1);
470  fgmres.Mult(B, U);
471  }
472  // 14. Recover the parallel grid function corresponding to U. This is the
473  // local finite element solution on each processor.
474  a->RecoverFEMSolution(U, b, u);
475 
476  if (exact_sol)
477  {
478  double err_r = -1.0;
479  double err_i = -1.0;
480 
481  switch (prob)
482  {
483  case 0:
484  err_r = u.real().ComputeL2Error(u0_r);
485  err_i = u.imag().ComputeL2Error(u0_i);
486  break;
487  case 1:
488  err_r = u.real().ComputeL2Error(u1_r);
489  err_i = u.imag().ComputeL2Error(u1_i);
490  break;
491  case 2:
492  err_r = u.real().ComputeL2Error(u2_r);
493  err_i = u.imag().ComputeL2Error(u2_i);
494  break;
495  default: break; // This should be unreachable
496  }
497 
498  if ( myid == 0 )
499  {
500  cout << endl;
501  cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
502  cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
503  cout << endl;
504  }
505  }
506 
507  // 15. Save the refined mesh and the solution in parallel. This output can be
508  // viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
509  {
510  ostringstream mesh_name, sol_r_name, sol_i_name;
511  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
512  sol_r_name << "sol_r." << setfill('0') << setw(6) << myid;
513  sol_i_name << "sol_i." << setfill('0') << setw(6) << myid;
514 
515  ofstream mesh_ofs(mesh_name.str().c_str());
516  mesh_ofs.precision(8);
517  pmesh->Print(mesh_ofs);
518 
519  ofstream sol_r_ofs(sol_r_name.str().c_str());
520  ofstream sol_i_ofs(sol_i_name.str().c_str());
521  sol_r_ofs.precision(8);
522  sol_i_ofs.precision(8);
523  u.real().Save(sol_r_ofs);
524  u.imag().Save(sol_i_ofs);
525  }
526 
527  // 16. Send the solution by socket to a GLVis server.
528  if (visualization)
529  {
530  char vishost[] = "localhost";
531  int visport = 19916;
532  socketstream sol_sock_r(vishost, visport);
533  socketstream sol_sock_i(vishost, visport);
534  sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
535  sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
536  sol_sock_r.precision(8);
537  sol_sock_i.precision(8);
538  sol_sock_r << "solution\n" << *pmesh << u.real()
539  << "window_title 'Solution: Real Part'" << flush;
540  sol_sock_i << "solution\n" << *pmesh << u.imag()
541  << "window_title 'Solution: Imaginary Part'" << flush;
542  }
543  if (visualization && exact_sol)
544  {
545  *u_exact -= u;
546 
547  char vishost[] = "localhost";
548  int visport = 19916;
549  socketstream sol_sock_r(vishost, visport);
550  socketstream sol_sock_i(vishost, visport);
551  sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
552  sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
553  sol_sock_r.precision(8);
554  sol_sock_i.precision(8);
555  sol_sock_r << "solution\n" << *pmesh << u_exact->real()
556  << "window_title 'Error: Real Part'" << flush;
557  sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
558  << "window_title 'Error: Imaginary Part'" << flush;
559  }
560  if (visualization)
561  {
562  ParGridFunction u_t(fespace);
563  u_t = u.real();
564  char vishost[] = "localhost";
565  int visport = 19916;
566  socketstream sol_sock(vishost, visport);
567  sol_sock << "parallel " << num_procs << " " << myid << "\n";
568  sol_sock.precision(8);
569  sol_sock << "solution\n" << *pmesh << u_t
570  << "window_title 'Harmonic Solution (t = 0.0 T)'"
571  << "pause\n" << flush;
572  if (myid == 0)
573  cout << "GLVis visualization paused."
574  << " Press space (in the GLVis window) to resume it.\n";
575  int num_frames = 32;
576  int i = 0;
577  while (sol_sock)
578  {
579  double t = (double)(i % num_frames) / num_frames;
580  ostringstream oss;
581  oss << "Harmonic Solution (t = " << t << " T)";
582 
583  add(cos( 2.0 * M_PI * t), u.real(),
584  sin(-2.0 * M_PI * t), u.imag(), u_t);
585  sol_sock << "parallel " << num_procs << " " << myid << "\n";
586  sol_sock << "solution\n" << *pmesh << u_t
587  << "window_title '" << oss.str() << "'" << flush;
588  i++;
589  }
590  }
591 
592  // 17. Free the used memory.
593  delete a;
594  delete u_exact;
595  delete pcOp;
596  delete fespace;
597  delete fec;
598  delete pmesh;
599 
600  return 0;
601 }
602 
603 bool check_for_inline_mesh(const char * mesh_file)
604 {
605  string file(mesh_file);
606  size_t p0 = file.find_last_of("/");
607  string s0 = file.substr((p0==string::npos)?0:(p0+1),7);
608  return s0 == "inline-";
609 }
610 
611 complex<double> u0_exact(const Vector &x)
612 {
613  int dim = x.Size();
614  complex<double> i(0.0, 1.0);
615  complex<double> alpha = (epsilon_ * omega_ - i * sigma_);
616  complex<double> kappa = std::sqrt(mu_ * omega_* alpha);
617  return std::exp(-i * kappa * x[dim - 1]);
618 }
619 
620 double u0_real_exact(const Vector &x)
621 {
622  return u0_exact(x).real();
623 }
624 
625 double u0_imag_exact(const Vector &x)
626 {
627  return u0_exact(x).imag();
628 }
629 
630 void u1_real_exact(const Vector &x, Vector &v)
631 {
632  int dim = x.Size();
633  v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
634 }
635 
636 void u1_imag_exact(const Vector &x, Vector &v)
637 {
638  int dim = x.Size();
639  v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
640 }
641 
642 void u2_real_exact(const Vector &x, Vector &v)
643 {
644  int dim = x.Size();
645  v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
646 }
647 
648 void u2_imag_exact(const Vector &x, Vector &v)
649 {
650  int dim = x.Size();
651  v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
652 }
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
int visport
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1743
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:1820
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
double u0_real_exact(const Vector &)
Definition: ex22p.cpp:620
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void u2_real_exact(const Vector &, Vector &)
Definition: ex22p.cpp:642
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
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.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1168
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
FGMRES method.
Definition: solvers.hpp:544
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
double u0_imag_exact(const Vector &)
Definition: ex22p.cpp:625
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
Abstract parallel finite element space.
Definition: pfespace.hpp:28
STL namespace.
(Q div u, div v) for RT elements
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:317
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
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. ...
char vishost[]
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
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void u2_imag_exact(const Vector &, Vector &)
Definition: ex22p.cpp:648
prob_type prob
Definition: ex25.cpp:154
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
A general vector function coefficient.
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
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
HYPRE_Int HYPRE_BigInt
int main(int argc, char *argv[])
Definition: ex22p.cpp:76
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
complex< double > u0_exact(const Vector &x)
Definition: ex22p.cpp:611
double a
Definition: lissajous.cpp:41
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
void u1_real_exact(const Vector &, Vector &)
Definition: ex22p.cpp:630
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
bool check_for_inline_mesh(const char *mesh_file)
Definition: ex22p.cpp:603
double freq
Definition: ex24.cpp:54
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Class for parallel bilinear form.
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 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
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4825
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Class for parallel grid function.
Definition: pgridfunc.hpp:32
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
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
void u1_imag_exact(const Vector &, Vector &)
Definition: ex22p.cpp:636
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).