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