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