MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex25p.cpp
Go to the documentation of this file.
1 // MFEM Example 25 - Parallel Version
2 //
3 // Compile with: make ex25p
4 //
5 // Sample runs: mpirun -np 4 ex25p -o 2 -f 1.0 -rs 1 -rp 1 -prob 0
6 // mpirun -np 4 ex25p -o 3 -f 10.0 -rs 1 -rp 1 -prob 1
7 // mpirun -np 4 ex25p -o 3 -f 5.0 -rs 3 -rp 1 -prob 2
8 // mpirun -np 4 ex25p -o 2 -f 1.0 -rs 1 -rp 1 -prob 3
9 // mpirun -np 4 ex25p -o 2 -f 1.0 -rs 2 -rp 2 -prob 0 -m ../data/beam-quad.mesh
10 // mpirun -np 4 ex25p -o 2 -f 8.0 -rs 2 -rp 2 -prob 4 -m ../data/inline-quad.mesh
11 // mpirun -np 4 ex25p -o 2 -f 2.0 -rs 1 -rp 1 -prob 4 -m ../data/inline-hex.mesh
12 //
13 // Device sample runs:
14 // mpirun -np 4 ex25p -o 1 -f 3.0 -rs 3 -rp 1 -prob 2 -pa -d cuda
15 // mpirun -np 4 ex25p -o 2 -f 1.0 -rs 1 -rp 1 -prob 3 -pa -d cuda
16 //
17 // Description: This example code solves a simple electromagnetic wave
18 // propagation problem corresponding to the second order
19 // indefinite Maxwell equation
20 //
21 // (1/mu) * curl curl E - \omega^2 * epsilon E = f
22 //
23 // with a Perfectly Matched Layer (PML).
24 //
25 // The example demonstrates discretization with Nedelec finite
26 // elements in 2D or 3D, as well as the use of complex-valued
27 // bilinear and linear forms. Several test problems are included,
28 // with prob = 0-3 having known exact solutions, see "On perfectly
29 // matched layers for discontinuous Petrov-Galerkin methods" by
30 // Vaziri Astaneh, Keith, Demkowicz, Comput Mech 63, 2019.
31 //
32 // We recommend viewing Example 22 before viewing this example.
33 
34 #include "mfem.hpp"
35 #include <fstream>
36 #include <iostream>
37 
38 #ifdef _WIN32
39 #define jn(n, x) _jn(n, x)
40 #define yn(n, x) _yn(n, x)
41 #endif
42 
43 using namespace std;
44 using namespace mfem;
45 
46 // Class for setting up a simple Cartesian PML region
47 class CartesianPML
48 {
49 private:
50  Mesh *mesh;
51 
52  int dim;
53 
54  // Length of the PML Region in each direction
55  Array2D<double> length;
56 
57  // Computational Domain Boundary
58  Array2D<double> comp_dom_bdr;
59 
60  // Domain Boundary
61  Array2D<double> dom_bdr;
62 
63  // Integer Array identifying elements in the PML
64  // 0: in the PML, 1: not in the PML
65  Array<int> elems;
66 
67  // Compute Domain and Computational Domain Boundaries
68  void SetBoundaries();
69 
70 public:
71  // Constructor
72  CartesianPML(Mesh *mesh_,Array2D<double> length_);
73 
74  // Return Computational Domain Boundary
75  Array2D<double> GetCompDomainBdr() {return comp_dom_bdr;}
76 
77  // Return Domain Boundary
78  Array2D<double> GetDomainBdr() {return dom_bdr;}
79 
80  // Return Markers list for elements
81  Array<int> * GetMarkedPMLElements() {return &elems;}
82 
83  // Mark elements in the PML region
84  void SetAttributes(ParMesh *pmesh);
85 
86  // PML complex stretching function
87  void StretchFunction(const Vector &x, vector<complex<double>> &dxs);
88 };
89 
90 // Class for returning the PML coefficients of the bilinear form
91 class PMLDiagMatrixCoefficient : public VectorCoefficient
92 {
93 private:
94  CartesianPML * pml = nullptr;
95  void (*Function)(const Vector &, CartesianPML * , Vector &);
96 public:
97  PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
98  Vector &),
99  CartesianPML * pml_)
100  : VectorCoefficient(dim), pml(pml_), Function(F)
101  {}
102 
103  using VectorCoefficient::Eval;
104 
105  virtual void Eval(Vector &K, ElementTransformation &T,
106  const IntegrationPoint &ip)
107  {
108  double x[3];
109  Vector transip(x, 3);
110  T.Transform(ip, transip);
111  K.SetSize(vdim);
112  (*Function)(transip, pml, K);
113  }
114 };
115 
116 void maxwell_solution(const Vector &x, vector<complex<double>> &Eval);
117 
118 void E_bdr_data_Re(const Vector &x, Vector &E);
119 void E_bdr_data_Im(const Vector &x, Vector &E);
120 
121 void E_exact_Re(const Vector &x, Vector &E);
122 void E_exact_Im(const Vector &x, Vector &E);
123 
124 void source(const Vector &x, Vector & f);
125 
126 // Functions for computing the necessary coefficients after PML stretching.
127 // J is the Jacobian matrix of the stretching function
128 void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, Vector & D);
129 void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, Vector & D);
130 void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, Vector & D);
131 
132 void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, Vector & D);
133 void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, Vector & D);
134 void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, Vector & D);
135 
138 
139 double mu = 1.0;
140 double epsilon = 1.0;
141 double omega;
142 int dim;
143 bool exact_known = false;
144 
146 {
147  beam, // Wave propagating in a beam-like domain
148  disc, // Point source propagating in the square-disc domain
149  lshape, // Point source propagating in the L-shape domain
150  fichera, // Point source propagating in the fichera domain
151  load_src // Approximated point source with PML all around
152 };
154 
155 int main(int argc, char *argv[])
156 {
157  // 1. Initialize MPI.
158  int num_procs, myid;
159  MPI_Init(&argc, &argv);
160  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
161  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
162 
163  // 2. Parse command-line options.
164  const char *mesh_file = nullptr;
165  int order = 1;
166  int ref_levels = 1;
167  int par_ref_levels = 2;
168  int iprob = 4;
169  double freq = 5.0;
170  bool herm_conv = true;
171  bool slu_solver = false;
172  bool mumps_solver = false;
173  bool visualization = 1;
174  bool pa = false;
175  const char *device_config = "cpu";
176 
177  OptionsParser args(argc, argv);
178  args.AddOption(&mesh_file, "-m", "--mesh",
179  "Mesh file to use.");
180  args.AddOption(&order, "-o", "--order",
181  "Finite element order (polynomial degree).");
182  args.AddOption(&iprob, "-prob", "--problem", "Problem case"
183  " 0: beam, 1: disc, 2: lshape, 3: fichera, 4: General");
184  args.AddOption(&ref_levels, "-rs", "--refinements-serial",
185  "Number of serial refinements");
186  args.AddOption(&par_ref_levels, "-rp", "--refinements-parallel",
187  "Number of parallel refinements");
188  args.AddOption(&mu, "-mu", "--permeability",
189  "Permeability of free space (or 1/(spring constant)).");
190  args.AddOption(&epsilon, "-eps", "--permittivity",
191  "Permittivity of free space (or mass constant).");
192  args.AddOption(&freq, "-f", "--frequency",
193  "Frequency (in Hz).");
194  args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
195  "--no-hermitian", "Use convention for Hermitian operators.");
196 #ifdef MFEM_USE_SUPERLU
197  args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
198  "--no-superlu", "Use the SuperLU Solver.");
199 #endif
200 #ifdef MFEM_USE_MUMPS
201  args.AddOption(&mumps_solver, "-mumps", "--mumps-solver", "-no-mumps",
202  "--no-mumps-solver", "Use the MUMPS Solver.");
203 #endif
204  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
205  "--no-visualization",
206  "Enable or disable GLVis visualization.");
207  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
208  "--no-partial-assembly", "Enable Partial Assembly.");
209  args.AddOption(&device_config, "-d", "--device",
210  "Device configuration string, see Device::Configure().");
211  args.Parse();
212  if (slu_solver && mumps_solver)
213  {
214  if (myid == 0)
215  cout << "WARNING: Both SuperLU and MUMPS have been selected,"
216  << " please choose either one." << endl
217  << " Defaulting to SuperLU." << endl;
218  mumps_solver = false;
219  }
220 
221  if (iprob > 4) { iprob = 4; }
222  prob = (prob_type)iprob;
223 
224  // 3. Enable hardware devices such as GPUs, and programming models such as
225  // CUDA, OCCA, RAJA and OpenMP based on command line options.
226  Device device(device_config);
227  if (myid == 0) { device.Print(); }
228 
229  // 4. Setup the (serial) mesh on all processors.
230  if (!mesh_file)
231  {
232  exact_known = true;
233  switch (prob)
234  {
235  case beam:
236  mesh_file = "../data/beam-hex.mesh";
237  break;
238  case disc:
239  mesh_file = "../data/square-disc.mesh";
240  break;
241  case lshape:
242  mesh_file = "../data/l-shape.mesh";
243  break;
244  case fichera:
245  mesh_file = "../data/fichera.mesh";
246  break;
247  default:
248  exact_known = false;
249  mesh_file = "../data/inline-quad.mesh";
250  break;
251  }
252  }
253 
254  if (!args.Good())
255  {
256  if (myid == 0)
257  {
258  args.PrintUsage(cout);
259  }
260  MPI_Finalize();
261  return 1;
262  }
263  if (myid == 0)
264  {
265  args.PrintOptions(cout);
266  }
267 
268  Mesh * mesh = new Mesh(mesh_file, 1, 1);
269  dim = mesh->Dimension();
270 
271  // Angular frequency
272  omega = 2.0 * M_PI * freq;
273 
274  // Setup PML length
275  Array2D<double> length(dim, 2); length = 0.0;
276 
277  // 5. Setup the Cartesian PML region.
278  switch (prob)
279  {
280  case disc:
281  length = 0.2;
282  break;
283  case lshape:
284  length(0, 0) = 0.1;
285  length(1, 0) = 0.1;
286  break;
287  case fichera:
288  length(0, 1) = 0.5;
289  length(1, 1) = 0.5;
290  length(2, 1) = 0.5;
291  break;
292  case beam:
293  length(0, 1) = 2.0;
294  break;
295  default:
296  length = 0.25;
297  break;
298  }
299  CartesianPML * pml = new CartesianPML(mesh,length);
300  comp_domain_bdr = pml->GetCompDomainBdr();
301  domain_bdr = pml->GetDomainBdr();
302 
303  // 6. Refine the serial mesh on all processors to increase the resolution.
304  for (int l = 0; l < ref_levels; l++)
305  {
306  mesh->UniformRefinement();
307  }
308 
309  // 7. Define a parallel mesh by a partitioning of the serial mesh.
310  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
311  delete mesh;
312  {
313  for (int l = 0; l < par_ref_levels; l++)
314  {
315  pmesh->UniformRefinement();
316  }
317  }
318 
319  // 7a. Reorient mesh in case of a tet mesh
320  pmesh->ReorientTetMesh();
321 
322  // 8. Set element attributes in order to distinguish elements in the PML
323  pml->SetAttributes(pmesh);
324 
325  // 9. Define a parallel finite element space on the parallel mesh. Here we
326  // use the Nedelec finite elements of the specified order.
327  FiniteElementCollection *fec = new ND_FECollection(order, dim);
328  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
329  HYPRE_BigInt size = fespace->GlobalTrueVSize();
330  if (myid == 0)
331  {
332  cout << "Number of finite element unknowns: " << size << endl;
333  }
334 
335  // 10. Determine the list of true (i.e. parallel conforming) essential
336  // boundary dofs. In this example, the boundary conditions are defined
337  // based on the specific mesh and the problem type.
338  Array<int> ess_tdof_list;
339  Array<int> ess_bdr;
340  if (pmesh->bdr_attributes.Size())
341  {
342  ess_bdr.SetSize(pmesh->bdr_attributes.Max());
343  ess_bdr = 1;
344  if (prob == lshape || prob == fichera)
345  {
346  ess_bdr = 0;
347  for (int j = 0; j < pmesh->GetNBE(); j++)
348  {
349  Vector center(dim);
350  int bdrgeom = pmesh->GetBdrElementBaseGeometry(j);
352  tr->Transform(Geometries.GetCenter(bdrgeom),center);
353  int k = pmesh->GetBdrAttribute(j);
354  switch (prob)
355  {
356  case lshape:
357  if (center[0] == 1.0 || center[0] == 0.5 || center[1] == 0.5)
358  {
359  ess_bdr[k - 1] = 1;
360  }
361  break;
362  case fichera:
363  if (center[0] == -1.0 || center[0] == 0.0 ||
364  center[1] == 0.0 || center[2] == 0.0)
365  {
366  ess_bdr[k - 1] = 1;
367  }
368  break;
369  default:
370  break;
371  }
372  }
373  }
374  }
375  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
376 
377  // 11. Setup Complex Operator convention
379  herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
380 
381  // 12. Set up the parallel linear form b(.) which corresponds to the
382  // right-hand side of the FEM linear system.
384  ParComplexLinearForm b(fespace, conv);
385  if (prob == load_src)
386  {
388  }
389  b.Vector::operator=(0.0);
390  b.Assemble();
391 
392  // 13. Define the solution vector x as a parallel complex finite element grid
393  // function corresponding to fespace.
394  ParComplexGridFunction x(fespace);
395  x = 0.0;
398  x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
399 
400  // 14. Set up the parallel sesquilinear form a(.,.)
401  //
402  // In Comp
403  // Domain: 1/mu (Curl E, Curl F) - omega^2 * epsilon (E,F)
404  //
405  // In PML: 1/mu (1/det(J) J^T J Curl E, Curl F)
406  // - omega^2 * epsilon (det(J) * (J^T J)^-1 * E, F)
407  //
408  // where J denotes the Jacobian Matrix of the PML Stretching function
409  Array<int> attr;
410  Array<int> attrPML;
411  if (pmesh->attributes.Size())
412  {
413  attr.SetSize(pmesh->attributes.Max());
414  attrPML.SetSize(pmesh->attributes.Max());
415  attr = 0; attr[0] = 1;
416  attrPML = 0;
417  if (pmesh->attributes.Max() > 1)
418  {
419  attrPML[1] = 1;
420  }
421  }
422 
423  ConstantCoefficient muinv(1.0/mu);
424  ConstantCoefficient omeg(-pow(omega, 2) * epsilon);
425  RestrictedCoefficient restr_muinv(muinv,attr);
426  RestrictedCoefficient restr_omeg(omeg,attr);
427 
428  // Integrators inside the computational domain (excluding the PML region)
429  ParSesquilinearForm a(fespace, conv);
430  a.AddDomainIntegrator(new CurlCurlIntegrator(restr_muinv),NULL);
431  a.AddDomainIntegrator(new VectorFEMassIntegrator(restr_omeg),NULL);
432 
433  int cdim = (dim == 2) ? 1 : dim;
434  PMLDiagMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, pml);
435  PMLDiagMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, pml);
436  ScalarVectorProductCoefficient c1_Re(muinv,pml_c1_Re);
437  ScalarVectorProductCoefficient c1_Im(muinv,pml_c1_Im);
438  VectorRestrictedCoefficient restr_c1_Re(c1_Re,attrPML);
439  VectorRestrictedCoefficient restr_c1_Im(c1_Im,attrPML);
440 
441  PMLDiagMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,pml);
442  PMLDiagMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,pml);
443  ScalarVectorProductCoefficient c2_Re(omeg,pml_c2_Re);
444  ScalarVectorProductCoefficient c2_Im(omeg,pml_c2_Im);
445  VectorRestrictedCoefficient restr_c2_Re(c2_Re,attrPML);
446  VectorRestrictedCoefficient restr_c2_Im(c2_Im,attrPML);
447 
448  // Integrators inside the PML region
449  a.AddDomainIntegrator(new CurlCurlIntegrator(restr_c1_Re),
450  new CurlCurlIntegrator(restr_c1_Im));
451  a.AddDomainIntegrator(new VectorFEMassIntegrator(restr_c2_Re),
452  new VectorFEMassIntegrator(restr_c2_Im));
453 
454  // 15. Assemble the parallel bilinear form and the corresponding linear
455  // system, applying any necessary transformations such as: parallel
456  // assembly, eliminating boundary conditions, applying conforming
457  // constraints for non-conforming AMR, etc.
458  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
459  a.Assemble();
460 
461  OperatorPtr Ah;
462  Vector B, X;
463  a.FormLinearSystem(ess_tdof_list, x, b, Ah, X, B);
464 
465  // 16. Solve using a direct or an iterative solver
466 #ifdef MFEM_USE_SUPERLU
467  if (!pa && slu_solver)
468  {
469  // Transform to monolithic HypreParMatrix
470  HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
471  SuperLURowLocMatrix SA(*A);
472  SuperLUSolver superlu(MPI_COMM_WORLD);
473  superlu.SetPrintStatistics(false);
474  superlu.SetSymmetricPattern(false);
476  superlu.SetOperator(SA);
477  superlu.Mult(B, X);
478  delete A;
479  }
480 #endif
481 #ifdef MFEM_USE_MUMPS
482  if (!pa && mumps_solver)
483  {
484  HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
485  MUMPSSolver mumps;
486  mumps.SetPrintLevel(0);
487  mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
488  mumps.SetOperator(*A);
489  mumps.Mult(B,X);
490  delete A;
491  }
492 #endif
493  // 16a. Set up the parallel Bilinear form a(.,.) for the preconditioner
494  //
495  // In Comp
496  // Domain: 1/mu (Curl E, Curl F) + omega^2 * epsilon (E,F)
497  //
498  // In PML: 1/mu (abs(1/det(J) J^T J) Curl E, Curl F)
499  // + omega^2 * epsilon (abs(det(J) * (J^T J)^-1) * E, F)
500  if (pa || (!slu_solver && !mumps_solver))
501  {
502  ConstantCoefficient absomeg(pow(omega, 2) * epsilon);
503  RestrictedCoefficient restr_absomeg(absomeg,attr);
504 
505  ParBilinearForm prec(fespace);
506  prec.AddDomainIntegrator(new CurlCurlIntegrator(restr_muinv));
507  prec.AddDomainIntegrator(new VectorFEMassIntegrator(restr_absomeg));
508 
509  PMLDiagMatrixCoefficient pml_c1_abs(cdim,detJ_inv_JT_J_abs, pml);
510  ScalarVectorProductCoefficient c1_abs(muinv,pml_c1_abs);
511  VectorRestrictedCoefficient restr_c1_abs(c1_abs,attrPML);
512 
513  PMLDiagMatrixCoefficient pml_c2_abs(dim, detJ_JT_J_inv_abs,pml);
514  ScalarVectorProductCoefficient c2_abs(absomeg,pml_c2_abs);
515  VectorRestrictedCoefficient restr_c2_abs(c2_abs,attrPML);
516 
517  prec.AddDomainIntegrator(new CurlCurlIntegrator(restr_c1_abs));
518  prec.AddDomainIntegrator(new VectorFEMassIntegrator(restr_c2_abs));
519 
520  if (pa) { prec.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
521  prec.Assemble();
522 
523  // 16b. Define and apply a parallel GMRES solver for AU=B with a block
524  // diagonal preconditioner based on hypre's AMS preconditioner.
525  Array<int> offsets(3);
526  offsets[0] = 0;
527  offsets[1] = fespace->GetTrueVSize();
528  offsets[2] = fespace->GetTrueVSize();
529  offsets.PartialSum();
530 
531  Operator *pc_r = nullptr;
532  Operator *pc_i = nullptr;
533  int s = (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0;
534  if (pa)
535  {
536  // Jacobi Smoother
537  OperatorJacobiSmoother *d00 = new OperatorJacobiSmoother(prec, ess_tdof_list);
538  ScaledOperator *d11 = new ScaledOperator(d00, s);
539  pc_r = d00;
540  pc_i = d11;
541  }
542  else
543  {
544  OperatorPtr PCOpAh;
545  prec.FormSystemMatrix(ess_tdof_list, PCOpAh);
546 
547  // Hypre AMS
548  HypreAMS *ams00 = new HypreAMS(*PCOpAh.As<HypreParMatrix>(), fespace);
549  ScaledOperator *ams11 = new ScaledOperator(ams00, s);
550  pc_r = ams00;
551  pc_i = ams11;
552  }
553 
554  BlockDiagonalPreconditioner BlockDP(offsets);
555  BlockDP.SetDiagonalBlock(0, pc_r);
556  BlockDP.SetDiagonalBlock(1, pc_i);
557 
558  GMRESSolver gmres(MPI_COMM_WORLD);
559  gmres.SetPrintLevel(1);
560  gmres.SetKDim(200);
561  gmres.SetMaxIter(pa ? 5000 : 2000);
562  gmres.SetRelTol(1e-5);
563  gmres.SetAbsTol(0.0);
564  gmres.SetOperator(*Ah);
565  gmres.SetPreconditioner(BlockDP);
566  gmres.Mult(B, X);
567  }
568 
569  // 17. Recover the parallel grid function corresponding to X. This is the
570  // local finite element solution on each processor.
571  a.RecoverFEMSolution(X, b, x);
572 
573  // If exact is known compute the error
574  if (exact_known)
575  {
578  int order_quad = max(2, 2 * order + 1);
579  const IntegrationRule *irs[Geometry::NumGeom];
580  for (int i = 0; i < Geometry::NumGeom; ++i)
581  {
582  irs[i] = &(IntRules.Get(i, order_quad));
583  }
584 
585  double L2Error_Re = x.real().ComputeL2Error(E_ex_Re, irs,
586  pml->GetMarkedPMLElements());
587  double L2Error_Im = x.imag().ComputeL2Error(E_ex_Im, irs,
588  pml->GetMarkedPMLElements());
589 
590  ParComplexGridFunction x_gf0(fespace);
591  x_gf0 = 0.0;
592  double norm_E_Re, norm_E_Im;
593  norm_E_Re = x_gf0.real().ComputeL2Error(E_ex_Re, irs,
594  pml->GetMarkedPMLElements());
595  norm_E_Im = x_gf0.imag().ComputeL2Error(E_ex_Im, irs,
596  pml->GetMarkedPMLElements());
597 
598  if (myid == 0)
599  {
600  cout << "\n Relative Error (Re part): || E_h - E || / ||E|| = "
601  << L2Error_Re / norm_E_Re
602  << "\n Relative Error (Im part): || E_h - E || / ||E|| = "
603  << L2Error_Im / norm_E_Im
604  << "\n Total Error: "
605  << sqrt(L2Error_Re*L2Error_Re + L2Error_Im*L2Error_Im) << "\n\n";
606  }
607  }
608 
609  // 18. Save the refined mesh and the solution in parallel. This output can be
610  // viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
611  {
612  ostringstream mesh_name, sol_r_name, sol_i_name;
613  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
614  sol_r_name << "ex25p-sol_r." << setfill('0') << setw(6) << myid;
615  sol_i_name << "ex25p-sol_i." << setfill('0') << setw(6) << myid;
616 
617  ofstream mesh_ofs(mesh_name.str().c_str());
618  mesh_ofs.precision(8);
619  pmesh->Print(mesh_ofs);
620 
621  ofstream sol_r_ofs(sol_r_name.str().c_str());
622  ofstream sol_i_ofs(sol_i_name.str().c_str());
623  sol_r_ofs.precision(8);
624  sol_i_ofs.precision(8);
625  x.real().Save(sol_r_ofs);
626  x.imag().Save(sol_i_ofs);
627  }
628 
629  // 19. Send the solution by socket to a GLVis server.
630  if (visualization)
631  {
632  // Define visualization keys for GLVis (see GLVis documentation)
633  string keys;
634  keys = (dim == 3) ? "keys macF\n" : keys = "keys amrRljcUUuu\n";
635  if (prob == beam && dim == 3) {keys = "keys macFFiYYYYYYYYYYYYYYYYYY\n";}
636  if (prob == beam && dim == 2) {keys = "keys amrRljcUUuuu\n"; }
637 
638  char vishost[] = "localhost";
639  int visport = 19916;
640 
641  {
642  socketstream sol_sock_re(vishost, visport);
643  sol_sock_re.precision(8);
644  sol_sock_re << "parallel " << num_procs << " " << myid << "\n"
645  << "solution\n" << *pmesh << x.real() << keys
646  << "window_title 'Solution real part'" << flush;
647  MPI_Barrier(MPI_COMM_WORLD); // try to prevent streams from mixing
648  }
649 
650  {
651  socketstream sol_sock_im(vishost, visport);
652  sol_sock_im.precision(8);
653  sol_sock_im << "parallel " << num_procs << " " << myid << "\n"
654  << "solution\n" << *pmesh << x.imag() << keys
655  << "window_title 'Solution imag part'" << flush;
656  MPI_Barrier(MPI_COMM_WORLD); // try to prevent streams from mixing
657  }
658 
659  {
660  ParGridFunction x_t(fespace);
661  x_t = x.real();
662 
663  socketstream sol_sock(vishost, visport);
664  sol_sock.precision(8);
665  sol_sock << "parallel " << num_procs << " " << myid << "\n"
666  << "solution\n" << *pmesh << x_t << keys << "autoscale off\n"
667  << "window_title 'Harmonic Solution (t = 0.0 T)'"
668  << "pause\n" << flush;
669 
670  if (myid == 0)
671  {
672  cout << "GLVis visualization paused."
673  << " Press space (in the GLVis window) to resume it.\n";
674  }
675 
676  int num_frames = 32;
677  int i = 0;
678  while (sol_sock)
679  {
680  double t = (double)(i % num_frames) / num_frames;
681  ostringstream oss;
682  oss << "Harmonic Solution (t = " << t << " T)";
683 
684  add(cos(2.0*M_PI*t), x.real(), sin(2.0*M_PI*t), x.imag(), x_t);
685  sol_sock << "parallel " << num_procs << " " << myid << "\n";
686  sol_sock << "solution\n" << *pmesh << x_t
687  << "window_title '" << oss.str() << "'" << flush;
688  i++;
689  }
690  }
691  }
692 
693  // 20. Free the used memory.
694  delete pml;
695  delete fespace;
696  delete fec;
697  delete pmesh;
698  MPI_Finalize();
699  return 0;
700 }
701 
702 void source(const Vector &x, Vector &f)
703 {
704  Vector center(dim);
705  double r = 0.0;
706  for (int i = 0; i < dim; ++i)
707  {
708  center(i) = 0.5 * (comp_domain_bdr(i, 0) + comp_domain_bdr(i, 1));
709  r += pow(x[i] - center[i], 2.);
710  }
711  double n = 5.0 * omega * sqrt(epsilon * mu) / M_PI;
712  double coeff = pow(n, 2) / M_PI;
713  double alpha = -pow(n, 2) * r;
714  f = 0.0;
715  f[0] = coeff * exp(alpha);
716 }
717 
718 void maxwell_solution(const Vector &x, vector<complex<double>> &E)
719 {
720  // Initialize
721  for (int i = 0; i < dim; ++i)
722  {
723  E[i] = 0.0;
724  }
725 
726  complex<double> zi = complex<double>(0., 1.);
727  double k = omega * sqrt(epsilon * mu);
728  switch (prob)
729  {
730  case disc:
731  case lshape:
732  case fichera:
733  {
734  Vector shift(dim);
735  shift = 0.0;
736  if (prob == fichera) { shift = 1.0; }
737  if (prob == disc) { shift = -0.5; }
738  if (prob == lshape) { shift = -1.0; }
739 
740  if (dim == 2)
741  {
742  double x0 = x(0) + shift(0);
743  double x1 = x(1) + shift(1);
744  double r = sqrt(x0 * x0 + x1 * x1);
745  double beta = k * r;
746 
747  // Bessel functions
748  complex<double> Ho, Ho_r, Ho_rr;
749  Ho = jn(0, beta) + zi * yn(0, beta);
750  Ho_r = -k * (jn(1, beta) + zi * yn(1, beta));
751  Ho_rr = -k * k * (1.0 / beta *
752  (jn(1, beta) + zi * yn(1, beta)) -
753  (jn(2, beta) + zi * yn(2, beta)));
754 
755  // First derivatives
756  double r_x = x0 / r;
757  double r_y = x1 / r;
758  double r_xy = -(r_x / r) * r_y;
759  double r_xx = (1.0 / r) * (1.0 - r_x * r_x);
760 
761  complex<double> val, val_xx, val_xy;
762  val = 0.25 * zi * Ho;
763  val_xx = 0.25 * zi * (r_xx * Ho_r + r_x * r_x * Ho_rr);
764  val_xy = 0.25 * zi * (r_xy * Ho_r + r_x * r_y * Ho_rr);
765  E[0] = zi / k * (k * k * val + val_xx);
766  E[1] = zi / k * val_xy;
767  }
768  else if (dim == 3)
769  {
770  double x0 = x(0) + shift(0);
771  double x1 = x(1) + shift(1);
772  double x2 = x(2) + shift(2);
773  double r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
774 
775  double r_x = x0 / r;
776  double r_y = x1 / r;
777  double r_z = x2 / r;
778  double r_xx = (1.0 / r) * (1.0 - r_x * r_x);
779  double r_yx = -(r_y / r) * r_x;
780  double r_zx = -(r_z / r) * r_x;
781 
782  complex<double> val, val_r, val_rr;
783  val = exp(zi * k * r) / r;
784  val_r = val / r * (zi * k * r - 1.0);
785  val_rr = val / (r * r) * (-k * k * r * r
786  - 2.0 * zi * k * r + 2.0);
787 
788  complex<double> val_xx, val_yx, val_zx;
789  val_xx = val_rr * r_x * r_x + val_r * r_xx;
790  val_yx = val_rr * r_x * r_y + val_r * r_yx;
791  val_zx = val_rr * r_x * r_z + val_r * r_zx;
792 
793  complex<double> alpha = zi * k / 4.0 / M_PI / k / k;
794  E[0] = alpha * (k * k * val + val_xx);
795  E[1] = alpha * val_yx;
796  E[2] = alpha * val_zx;
797  }
798  break;
799  }
800  case beam:
801  {
802  // T_10 mode
803  if (dim == 3)
804  {
805  double k10 = sqrt(k * k - M_PI * M_PI);
806  E[1] = -zi * k / M_PI * sin(M_PI*x(2))*exp(zi * k10 * x(0));
807  }
808  else if (dim == 2)
809  {
810  E[1] = -zi * k / M_PI * exp(zi * k * x(0));
811  }
812  break;
813  }
814  default:
815  break;
816  }
817 }
818 
819 void E_exact_Re(const Vector &x, Vector &E)
820 {
821  vector<complex<double>> Eval(E.Size());
822  maxwell_solution(x, Eval);
823  for (int i = 0; i < dim; ++i)
824  {
825  E[i] = Eval[i].real();
826  }
827 }
828 
829 void E_exact_Im(const Vector &x, Vector &E)
830 {
831  vector<complex<double>> Eval(E.Size());
832  maxwell_solution(x, Eval);
833  for (int i = 0; i < dim; ++i)
834  {
835  E[i] = Eval[i].imag();
836  }
837 }
838 
839 void E_bdr_data_Re(const Vector &x, Vector &E)
840 {
841  E = 0.0;
842  bool in_pml = false;
843 
844  for (int i = 0; i < dim; ++i)
845  {
846  // check if in PML
847  if (x(i) - comp_domain_bdr(i, 0) < 0.0 ||
848  x(i) - comp_domain_bdr(i, 1) > 0.0)
849  {
850  in_pml = true;
851  break;
852  }
853  }
854  if (!in_pml)
855  {
856  vector<complex<double>> Eval(E.Size());
857  maxwell_solution(x, Eval);
858  for (int i = 0; i < dim; ++i)
859  {
860  E[i] = Eval[i].real();
861  }
862  }
863 }
864 
865 // Define bdr_data solution
866 void E_bdr_data_Im(const Vector &x, Vector &E)
867 {
868  E = 0.0;
869  bool in_pml = false;
870 
871  for (int i = 0; i < dim; ++i)
872  {
873  // check if in PML
874  if (x(i) - comp_domain_bdr(i, 0) < 0.0 ||
875  x(i) - comp_domain_bdr(i, 1) > 0.0)
876  {
877  in_pml = true;
878  break;
879  }
880  }
881  if (!in_pml)
882  {
883  vector<complex<double>> Eval(E.Size());
884  maxwell_solution(x, Eval);
885  for (int i = 0; i < dim; ++i)
886  {
887  E[i] = Eval[i].imag();
888  }
889  }
890 }
891 
892 void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, Vector & D)
893 {
894  vector<complex<double>> dxs(dim);
895  complex<double> det(1.0, 0.0);
896  pml->StretchFunction(x, dxs);
897 
898  for (int i = 0; i < dim; ++i)
899  {
900  det *= dxs[i];
901  }
902 
903  for (int i = 0; i < dim; ++i)
904  {
905  D(i) = (det / pow(dxs[i], 2)).real();
906  }
907 }
908 
909 void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, Vector & D)
910 {
911  vector<complex<double>> dxs(dim);
912  complex<double> det = 1.0;
913  pml->StretchFunction(x, dxs);
914 
915  for (int i = 0; i < dim; ++i)
916  {
917  det *= dxs[i];
918  }
919 
920  for (int i = 0; i < dim; ++i)
921  {
922  D(i) = (det / pow(dxs[i], 2)).imag();
923  }
924 }
925 
926 void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, Vector & D)
927 {
928  vector<complex<double>> dxs(dim);
929  complex<double> det = 1.0;
930  pml->StretchFunction(x, dxs);
931 
932  for (int i = 0; i < dim; ++i)
933  {
934  det *= dxs[i];
935  }
936 
937  for (int i = 0; i < dim; ++i)
938  {
939  D(i) = abs(det / pow(dxs[i], 2));
940  }
941 }
942 
943 void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, Vector & D)
944 {
945  vector<complex<double>> dxs(dim);
946  complex<double> det(1.0, 0.0);
947  pml->StretchFunction(x, dxs);
948 
949  for (int i = 0; i < dim; ++i)
950  {
951  det *= dxs[i];
952  }
953 
954  // in the 2D case the coefficient is scalar 1/det(J)
955  if (dim == 2)
956  {
957  D = (1.0 / det).real();
958  }
959  else
960  {
961  for (int i = 0; i < dim; ++i)
962  {
963  D(i) = (pow(dxs[i], 2) / det).real();
964  }
965  }
966 }
967 
968 void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, Vector & D)
969 {
970  vector<complex<double>> dxs(dim);
971  complex<double> det = 1.0;
972  pml->StretchFunction(x, dxs);
973 
974  for (int i = 0; i < dim; ++i)
975  {
976  det *= dxs[i];
977  }
978 
979  if (dim == 2)
980  {
981  D = (1.0 / det).imag();
982  }
983  else
984  {
985  for (int i = 0; i < dim; ++i)
986  {
987  D(i) = (pow(dxs[i], 2) / det).imag();
988  }
989  }
990 }
991 
992 void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, Vector & D)
993 {
994  vector<complex<double>> dxs(dim);
995  complex<double> det = 1.0;
996  pml->StretchFunction(x, dxs);
997 
998  for (int i = 0; i < dim; ++i)
999  {
1000  det *= dxs[i];
1001  }
1002 
1003  if (dim == 2)
1004  {
1005  D = abs(1.0 / det);
1006  }
1007  else
1008  {
1009  for (int i = 0; i < dim; ++i)
1010  {
1011  D(i) = abs(pow(dxs[i], 2) / det);
1012  }
1013  }
1014 }
1015 
1016 CartesianPML::CartesianPML(Mesh *mesh_, Array2D<double> length_)
1017  : mesh(mesh_), length(length_)
1018 {
1019  dim = mesh->Dimension();
1020  SetBoundaries();
1021 }
1022 
1023 void CartesianPML::SetBoundaries()
1024 {
1025  comp_dom_bdr.SetSize(dim, 2);
1026  dom_bdr.SetSize(dim, 2);
1027  Vector pmin, pmax;
1028  mesh->GetBoundingBox(pmin, pmax);
1029  for (int i = 0; i < dim; i++)
1030  {
1031  dom_bdr(i, 0) = pmin(i);
1032  dom_bdr(i, 1) = pmax(i);
1033  comp_dom_bdr(i, 0) = dom_bdr(i, 0) + length(i, 0);
1034  comp_dom_bdr(i, 1) = dom_bdr(i, 1) - length(i, 1);
1035  }
1036 }
1037 
1038 void CartesianPML::SetAttributes(ParMesh *pmesh)
1039 {
1040  // Initialize bdr attributes
1041  for (int i = 0; i < pmesh->GetNBE(); ++i)
1042  {
1043  pmesh->GetBdrElement(i)->SetAttribute(i+1);
1044  }
1045 
1046  int nrelem = pmesh->GetNE();
1047 
1048  // Initialize list with 1
1049  elems.SetSize(nrelem);
1050 
1051  // Loop through the elements and identify which of them are in the PML
1052  for (int i = 0; i < nrelem; ++i)
1053  {
1054  elems[i] = 1;
1055  bool in_pml = false;
1056  Element *el = pmesh->GetElement(i);
1057  Array<int> vertices;
1058 
1059  // Initialize attribute
1060  el->SetAttribute(1);
1061  el->GetVertices(vertices);
1062  int nrvert = vertices.Size();
1063 
1064  // Check if any vertex is in the PML
1065  for (int iv = 0; iv < nrvert; ++iv)
1066  {
1067  int vert_idx = vertices[iv];
1068  double *coords = pmesh->GetVertex(vert_idx);
1069  for (int comp = 0; comp < dim; ++comp)
1070  {
1071  if (coords[comp] > comp_dom_bdr(comp, 1) ||
1072  coords[comp] < comp_dom_bdr(comp, 0))
1073  {
1074  in_pml = true;
1075  break;
1076  }
1077  }
1078  }
1079  if (in_pml)
1080  {
1081  elems[i] = 0;
1082  el->SetAttribute(2);
1083  }
1084  }
1085  pmesh->SetAttributes();
1086 }
1087 
1088 void CartesianPML::StretchFunction(const Vector &x,
1089  vector<complex<double>> &dxs)
1090 {
1091  complex<double> zi = complex<double>(0., 1.);
1092 
1093  double n = 2.0;
1094  double c = 5.0;
1095  double coeff;
1096  double k = omega * sqrt(epsilon * mu);
1097 
1098  // Stretch in each direction independently
1099  for (int i = 0; i < dim; ++i)
1100  {
1101  dxs[i] = 1.0;
1102  if (x(i) >= comp_domain_bdr(i, 1))
1103  {
1104  coeff = n * c / k / pow(length(i, 1), n);
1105  dxs[i] = 1.0 + zi * coeff *
1106  abs(pow(x(i) - comp_domain_bdr(i, 1), n - 1.0));
1107  }
1108  if (x(i) <= comp_domain_bdr(i, 0))
1109  {
1110  coeff = n * c / k / pow(length(i, 0), n);
1111  dxs[i] = 1.0 + zi * coeff *
1112  abs(pow(x(i) - comp_domain_bdr(i, 0), n - 1.0));
1113  }
1114  }
1115 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:790
Definition: ex25.cpp:148
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:99
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: superlu.cpp:474
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1203
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1540
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:917
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:920
void SetOperator(const Operator &op)
Set the Operator and perform factorization.
Definition: mumps.cpp:30
Base class for vector Coefficients that optionally depend on time and space.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
double epsilon
Definition: ex25.cpp:140
void E_exact_Re(const Vector &x, Vector &E)
Definition: ex25.cpp:739
virtual void ReorientTetMesh()
See the remarks for the serial version in mesh.hpp.
Definition: pmesh.cpp:2798
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:849
Integrator for (curl u, curl v) for Nedelec elements.
void SetPrintLevel(int print_lvl)
Set the error print level for MUMPS.
Definition: mumps.cpp:247
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
void maxwell_solution(const Vector &x, vector< complex< double >> &Eval)
Definition: ex25.cpp:638
void AddDomainIntegrator(LinearFormIntegrator *lfi_real, LinearFormIntegrator *lfi_imag)
Adds new Domain Integrator.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
void detJ_JT_J_inv_abs(const Vector &x, CartesianPML *pml, Vector &D)
Definition: ex25.cpp:846
ParGridFunction & imag()
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:841
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.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void detJ_inv_JT_J_abs(const Vector &x, CartesianPML *pml, Vector &D)
Definition: ex25.cpp:912
void E_bdr_data_Re(const Vector &x, Vector &E)
Definition: ex25.cpp:759
virtual void SetAttributes()
Definition: pmesh.cpp:1574
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:70
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:291
void SetSymmetricPattern(bool sym)
Definition: superlu.cpp:423
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:279
ElementTransformation * GetBdrElementTransformation(int i)
Returns the transformation defining the i-th boundary element.
Definition: mesh.cpp:428
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
void Mult(const Vector &x, Vector &y) const
Solve y = Op^{-1} x.
Definition: mumps.cpp:197
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: superlu.cpp:580
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:884
void SetMatrixSymType(MatType mtype)
Set the matrix type.
Definition: mumps.cpp:252
Derived vector coefficient that has the value of the parent vector where it is active and is zero oth...
bool exact_known
Definition: ex25.cpp:143
void source(const Vector &x, Vector &f)
Definition: ex25.cpp:622
Geometry Geometries
Definition: fe.cpp:13507
void E_exact_Im(const Vector &x, Vector &E)
Definition: ex25.cpp:749
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
Array2D< double > domain_bdr
Definition: ex25.cpp:137
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. ...
constexpr char vishost[]
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
Definition: ex25.cpp:149
void E_bdr_data_Im(const Vector &x, Vector &E)
Definition: ex25.cpp:786
void SetColumnPermutation(superlu::ColPerm col_perm)
Definition: superlu.cpp:343
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
void detJ_inv_JT_J_Im(const Vector &x, CartesianPML *pml, Vector &D)
Definition: ex25.cpp:888
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
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:361
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4382
void SetPrintStatistics(bool print_stat)
Definition: superlu.cpp:325
const Element * GetElement(int i) const
Definition: mesh.hpp:942
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
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
A general vector function coefficient.
double omega
Definition: ex25.cpp:141
void SetAbsTol(double atol)
Definition: solvers.hpp:99
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:204
Array2D< double > comp_domain_bdr
Definition: ex25.cpp:136
void SetRelTol(double rtol)
Definition: solvers.hpp:98
Scaled Operator B: x -&gt; a A(x).
Definition: operator.hpp:692
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
Vector coefficient defined as a product of scalar and vector coefficients.
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
HYPRE_Int HYPRE_BigInt
GMRES method.
Definition: solvers.hpp:348
void AddDomainIntegrator(BilinearFormIntegrator *bfi_real, BilinearFormIntegrator *bfi_imag)
Adds new Domain Integrator.
void detJ_JT_J_inv_Re(const Vector &x, CartesianPML *pml, Vector &D)
Definition: ex25.cpp:812
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 Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff, VectorCoefficient &imag_coeff, Array< int > &attr)
Class for integration point with weight.
Definition: intrules.hpp:25
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:273
double mu
Definition: ex25.cpp:139
int dim
Definition: ex24.cpp:53
void detJ_JT_J_inv_Im(const Vector &x, CartesianPML *pml, Vector &D)
Definition: ex25.cpp:829
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
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 SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:58
Class for parallel bilinear form.
Geometry::Type GetBdrElementBaseGeometry(int i) const
Definition: mesh.hpp:974
RefCoord t[3]
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:266
const double alpha
Definition: ex15.cpp:369
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:416
Vector data type.
Definition: vector.hpp:60
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Definition: ex25.cpp:147
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
ParGridFunction & real()
RefCoord s[3]
Class for parallel grid function.
Definition: pgridfunc.hpp:32
prob_type
Definition: ex25.cpp:145
void detJ_inv_JT_J_Re(const Vector &x, CartesianPML *pml, Vector &D)
Definition: ex25.cpp:863
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:277
Class for parallel meshes.
Definition: pmesh.hpp:32
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:377
Abstract data type element.
Definition: element.hpp:28
int main()
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:202
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:946
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
MUMPS: A Parallel Sparse Direct Solver.
Definition: mumps.hpp:34
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150