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