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