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