MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex25.cpp
Go to the documentation of this file.
1// MFEM Example 25
2//
3// Compile with: make ex25
4//
5// Sample runs: ex25 -o 2 -f 1.0 -ref 2 -prob 0
6// ex25 -o 3 -f 10.0 -ref 2 -prob 1
7// ex25 -o 2 -f 5.0 -ref 4 -prob 2
8// ex25 -o 2 -f 1.0 -ref 2 -prob 3
9// ex25 -o 2 -f 1.0 -ref 2 -prob 0 -m ../data/beam-quad.mesh
10// ex25 -o 2 -f 8.0 -ref 3 -prob 4 -m ../data/inline-quad.mesh
11// ex25 -o 2 -f 2.0 -ref 1 -prob 4 -m ../data/inline-hex.mesh
12//
13// Device sample runs:
14// ex25 -o 2 -f 8.0 -ref 3 -prob 4 -m ../data/inline-quad.mesh -pa -d cuda
15// ex25 -o 2 -f 2.0 -ref 1 -prob 4 -m ../data/inline-hex.mesh -pa -d cuda
16//
17// Description: This example code solves a simple electromagnetic wave
18// propagation problem corresponding to the second order
19// indefinite Maxwell equation
20//
21// (1/mu) * curl curl E - \omega^2 * epsilon E = f
22//
23// with a Perfectly Matched Layer (PML).
24//
25// The example demonstrates discretization with Nedelec finite
26// elements in 2D or 3D, as well as the use of complex-valued
27// bilinear and linear forms. Several test problems are included,
28// with prob = 0-3 having known exact solutions, see "On perfectly
29// matched layers for discontinuous Petrov-Galerkin methods" by
30// Vaziri Astaneh, Keith, Demkowicz, Comput Mech 63, 2019.
31//
32// We recommend viewing Example 22 before viewing this example.
33
34#include "mfem.hpp"
35#include <memory>
36#include <fstream>
37#include <iostream>
38
39#ifdef _WIN32
40#define jn(n, x) _jn(n, x)
41#define yn(n, x) _yn(n, x)
42#endif
43
44using namespace std;
45using namespace mfem;
46
47// Class for setting up a simple Cartesian PML region
48class PML
49{
50private:
51 Mesh *mesh;
52
53 int dim;
54
55 // Length of the PML Region in each direction
56 Array2D<real_t> length;
57
58 // Computational Domain Boundary
59 Array2D<real_t> comp_dom_bdr;
60
61 // Domain Boundary
62 Array2D<real_t> dom_bdr;
63
64 // Integer Array identifying elements in the PML
65 // 0: in the PML, 1: not in the PML
66 Array<int> elems;
67
68 // Compute Domain and Computational Domain Boundaries
69 void SetBoundaries();
70
71public:
72 // Constructor
73 PML(Mesh *mesh_,Array2D<real_t> length_);
74
75 // Return Computational Domain Boundary
76 Array2D<real_t> GetCompDomainBdr() {return comp_dom_bdr;}
77
78 // Return Domain Boundary
79 Array2D<real_t> GetDomainBdr() {return dom_bdr;}
80
81 // Return Markers list for elements
82 Array<int> * GetMarkedPMLElements() {return &elems;}
83
84 // Mark elements in the PML region
85 void SetAttributes(Mesh *mesh_);
86
87 // PML complex stretching function
88 void StretchFunction(const Vector &x, vector<complex<real_t>> &dxs);
89};
90
91// Class for returning the PML coefficients of the bilinear form
92class PMLDiagMatrixCoefficient : public VectorCoefficient
93{
94private:
95 PML * pml = nullptr;
96 void (*Function)(const Vector &, PML *, Vector &);
97public:
98 PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, PML *,
99 Vector &),
100 PML * pml_)
101 : VectorCoefficient(dim), pml(pml_), Function(F)
102 {}
103
105
106 virtual void Eval(Vector &K, ElementTransformation &T,
107 const IntegrationPoint &ip)
108 {
109 real_t x[3];
110 Vector transip(x, 3);
111 T.Transform(ip, transip);
112 K.SetSize(vdim);
113 (*Function)(transip, pml, K);
114 }
115};
116
117void maxwell_solution(const Vector &x, vector<complex<real_t>> &Eval);
118
119void E_bdr_data_Re(const Vector &x, Vector &E);
120void E_bdr_data_Im(const Vector &x, Vector &E);
121
122void E_exact_Re(const Vector &x, Vector &E);
123void E_exact_Im(const Vector &x, Vector &E);
124
125void source(const Vector &x, Vector & f);
126
127// Functions for computing the necessary coefficients after PML stretching.
128// J is the Jacobian matrix of the stretching function
129void detJ_JT_J_inv_Re(const Vector &x, PML * pml, Vector &D);
130void detJ_JT_J_inv_Im(const Vector &x, PML * pml, Vector &D);
131void detJ_JT_J_inv_abs(const Vector &x, PML * pml, Vector &D);
132
133void detJ_inv_JT_J_Re(const Vector &x, PML * pml, Vector &D);
134void detJ_inv_JT_J_Im(const Vector &x, PML * pml, Vector &D);
135void detJ_inv_JT_J_abs(const Vector &x, PML * pml, Vector &D);
136
139
140real_t mu = 1.0;
143int dim;
144bool exact_known = false;
145
146template <typename T> T pow2(const T &x) { return x*x; }
147
149{
150 beam, // Wave propagating in a beam-like domain
151 disc, // Point source propagating in the square-disc domain
152 lshape, // Point source propagating in the L-shape domain
153 fichera, // Point source propagating in the fichera domain
154 load_src // Approximated point source with PML all around
157
158int main(int argc, char *argv[])
159{
160 // 1. Parse command-line options.
161 const char *mesh_file = nullptr;
162 int order = 1;
163 int ref_levels = 3;
164 int iprob = 4;
165 real_t freq = 5.0;
166 bool herm_conv = true;
167 bool umf_solver = false;
168 bool visualization = 1;
169 bool pa = false;
170 const char *device_config = "cpu";
171
172 OptionsParser args(argc, argv);
174 "Mesh file to use.");
176 "Finite element order (polynomial degree).");
177 args.AddOption(&iprob, "-prob", "--problem", "Problem case"
178 " 0: beam, 1: disc, 2: lshape, 3: fichera, 4: General");
180 "Number of refinements");
182 "Permeability of free space (or 1/(spring constant)).");
184 "Permittivity of free space (or mass constant).");
186 "Frequency (in Hz).");
188 "--no-hermitian", "Use convention for Hermitian operators.");
189#ifdef MFEM_USE_SUITESPARSE
191 "--no-umfpack", "Use the UMFPack Solver.");
192#endif
194 "--no-visualization",
195 "Enable or disable GLVis visualization.");
197 "--no-partial-assembly", "Enable Partial Assembly.");
199 "Device configuration string, see Device::Configure().");
200 args.Parse();
201
202 if (iprob > 4) { iprob = 4; }
203 prob = (prob_type)iprob;
204
205 // 2. Enable hardware devices such as GPUs, and programming models such as
206 // CUDA, OCCA, RAJA and OpenMP based on command line options.
207 Device device(device_config);
208 device.Print();
209
210 // 3. Setup the mesh
211 if (!mesh_file)
212 {
213 exact_known = true;
214 switch (prob)
215 {
216 case beam:
217 mesh_file = "../data/beam-hex.mesh";
218 break;
219 case disc:
220 mesh_file = "../data/square-disc.mesh";
221 break;
222 case lshape:
223 mesh_file = "../data/l-shape.mesh";
224 break;
225 case fichera:
226 mesh_file = "../data/fichera.mesh";
227 break;
228 default:
229 exact_known = false;
231 break;
232 }
233 }
234
235 if (!args.Good())
236 {
237 args.PrintUsage(cout);
238 return 1;
239 }
240 args.PrintOptions(cout);
241
242 Mesh * mesh = new Mesh(mesh_file, 1, 1);
243 dim = mesh->Dimension();
244
245 // Angular frequency
246 omega = real_t(2.0 * M_PI) * freq;
247
248 // Setup PML length
249 Array2D<real_t> length(dim, 2); length = 0.0;
250
251 // 4. Setup the Cartesian PML region.
252 switch (prob)
253 {
254 case disc:
255 length = 0.2;
256 break;
257 case lshape:
258 length(0, 0) = 0.1;
259 length(1, 0) = 0.1;
260 break;
261 case fichera:
262 length(0, 1) = 0.5;
263 length(1, 1) = 0.5;
264 length(2, 1) = 0.5;
265 break;
266 case beam:
267 length(0, 1) = 2.0;
268 break;
269 default:
270 length = 0.25;
271 break;
272 }
273 PML * pml = new PML(mesh,length);
274 comp_domain_bdr = pml->GetCompDomainBdr();
275 domain_bdr = pml->GetDomainBdr();
276
277 // 5. Refine the mesh to increase the resolution.
278 for (int l = 0; l < ref_levels; l++)
279 {
280 mesh->UniformRefinement();
281 }
282
283 // 6. Set element attributes in order to distinguish elements in the
284 // PML region
285 pml->SetAttributes(mesh);
286
287 // 7. Define a finite element space on the mesh. Here we use the Nedelec
288 // finite elements of the specified order.
290 FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
291 int size = fespace->GetTrueVSize();
292
293 cout << "Number of finite element unknowns: " << size << endl;
294
295 // 8. Determine the list of true essential boundary dofs. In this example,
296 // the boundary conditions are defined based on the specific mesh and the
297 // problem type.
298 Array<int> ess_tdof_list;
299 Array<int> ess_bdr;
300 if (mesh->bdr_attributes.Size())
301 {
302 ess_bdr.SetSize(mesh->bdr_attributes.Max());
303 ess_bdr = 1;
304 if (prob == lshape || prob == fichera)
305 {
306 ess_bdr = 0;
307 for (int j = 0; j < mesh->GetNBE(); j++)
308 {
309 Vector center(dim);
310 int bdrgeom = mesh->GetBdrElementGeometry(j);
312 tr->Transform(Geometries.GetCenter(bdrgeom),center);
313 int k = mesh->GetBdrAttribute(j);
314 switch (prob)
315 {
316 case lshape:
317 if (center[0] == 1_r || center[0] == 0.5_r ||
318 center[1] == 0.5_r)
319 {
320 ess_bdr[k - 1] = 1;
321 }
322 break;
323 case fichera:
324 if (center[0] == -1_r || center[0] == 0_r ||
325 center[1] == 0_r || center[2] == 0_r)
326 {
327 ess_bdr[k - 1] = 1;
328 }
329 break;
330 default:
331 break;
332 }
333 }
334 }
335 }
336 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
337
338 // 9. Setup Complex Operator convention
341
342 // 10. Set up the linear form b(.) which corresponds to the right-hand side of
343 // the FEM linear system.
345 ComplexLinearForm b(fespace, conv);
347 {
349 }
350 b.Vector::operator=(0.0);
351 b.Assemble();
352
353 // 11. Define the solution vector x as a complex finite element grid function
354 // corresponding to fespace.
355 ComplexGridFunction x(fespace);
356 x = 0.0;
359 x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
360
361 // 12. Set up the sesquilinear form a(.,.)
362 //
363 // In Comp
364 // Domain: 1/mu (Curl E, Curl F) - omega^2 * epsilon (E,F)
365 //
366 // In PML: 1/mu (1/det(J) J^T J Curl E, Curl F)
367 // - omega^2 * epsilon (det(J) * (J^T J)^-1 * E, F)
368 //
369 // where J denotes the Jacobian Matrix of the PML Stretching function
370 Array<int> attr;
371 Array<int> attrPML;
372 if (mesh->attributes.Size())
373 {
374 attr.SetSize(mesh->attributes.Max());
375 attrPML.SetSize(mesh->attributes.Max());
376 attr = 0; attr[0] = 1;
377 attrPML = 0;
378 if (mesh->attributes.Max() > 1)
379 {
380 attrPML[1] = 1;
381 }
382 }
383
384 ConstantCoefficient muinv(1_r / mu);
386 RestrictedCoefficient restr_muinv(muinv,attr);
387 RestrictedCoefficient restr_omeg(omeg,attr);
388
389 // Integrators inside the computational domain (excluding the PML region)
390 SesquilinearForm a(fespace, conv);
393
394 int cdim = (dim == 2) ? 1 : dim;
395 PMLDiagMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, pml);
396 PMLDiagMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, pml);
397 ScalarVectorProductCoefficient c1_Re(muinv,pml_c1_Re);
398 ScalarVectorProductCoefficient c1_Im(muinv,pml_c1_Im);
399 VectorRestrictedCoefficient restr_c1_Re(c1_Re,attrPML);
400 VectorRestrictedCoefficient restr_c1_Im(c1_Im,attrPML);
401
402 PMLDiagMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,pml);
403 PMLDiagMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,pml);
404 ScalarVectorProductCoefficient c2_Re(omeg,pml_c2_Re);
405 ScalarVectorProductCoefficient c2_Im(omeg,pml_c2_Im);
406 VectorRestrictedCoefficient restr_c2_Re(c2_Re,attrPML);
407 VectorRestrictedCoefficient restr_c2_Im(c2_Im,attrPML);
408
409 // Integrators inside the PML region
411 new CurlCurlIntegrator(restr_c1_Im));
413 new VectorFEMassIntegrator(restr_c2_Im));
414
415 // 13. Assemble the bilinear form and the corresponding linear system,
416 // applying any necessary transformations such as: assembly, eliminating
417 // boundary conditions, applying conforming constraints for
418 // non-conforming AMR, etc.
419 if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
420 a.Assemble(0);
421
422 OperatorPtr A;
423 Vector B, X;
424 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
425
426 // 14. Solve using a direct or an iterative solver
427#ifdef MFEM_USE_SUITESPARSE
428 if (!pa && umf_solver)
429 {
431 csolver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
432 csolver.SetPrintLevel(1);
433 csolver.Mult(B, X);
434 }
435#endif
436 // 14a. Set up the Bilinear form a(.,.) for the preconditioner
437 //
438 // In Comp
439 // Domain: 1/mu (Curl E, Curl F) + omega^2 * epsilon (E,F)
440 //
441 // In PML: 1/mu (abs(1/det(J) J^T J) Curl E, Curl F)
442 // + omega^2 * epsilon (abs(det(J) * (J^T J)^-1) * E, F)
443 if (pa || !umf_solver)
444 {
446 RestrictedCoefficient restr_absomeg(absomeg,attr);
447
448 BilinearForm prec(fespace);
451
452 PMLDiagMatrixCoefficient pml_c1_abs(cdim,detJ_inv_JT_J_abs, pml);
453 ScalarVectorProductCoefficient c1_abs(muinv,pml_c1_abs);
454 VectorRestrictedCoefficient restr_c1_abs(c1_abs,attrPML);
455
456 PMLDiagMatrixCoefficient pml_c2_abs(dim, detJ_JT_J_inv_abs,pml);
457 ScalarVectorProductCoefficient c2_abs(absomeg,pml_c2_abs);
458 VectorRestrictedCoefficient restr_c2_abs(c2_abs,attrPML);
459
462
463 if (pa) { prec.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
464 prec.Assemble();
465
466 // 14b. Define and apply a GMRES solver for AU=B with a block diagonal
467 // preconditioner based on the Gauss-Seidel or Jacobi sparse smoother.
468 Array<int> offsets(3);
469 offsets[0] = 0;
470 offsets[1] = fespace->GetTrueVSize();
471 offsets[2] = fespace->GetTrueVSize();
472 offsets.PartialSum();
473
474 std::unique_ptr<Operator> pc_r;
475 std::unique_ptr<Operator> pc_i;
476 real_t s = (conv == ComplexOperator::HERMITIAN) ? -1_r : 1_r;
477 if (pa)
478 {
479 // Jacobi Smoother
480 pc_r.reset(new OperatorJacobiSmoother(prec, ess_tdof_list));
481 pc_i.reset(new ScaledOperator(pc_r.get(), s));
482 }
483 else
484 {
485 OperatorPtr PCOpAh;
487 prec.FormSystemMatrix(ess_tdof_list, PCOpAh);
488
489 // Gauss-Seidel Smoother
490 pc_r.reset(new GSSmoother(*PCOpAh.As<SparseMatrix>()));
491 pc_i.reset(new ScaledOperator(pc_r.get(), s));
492 }
493
494 BlockDiagonalPreconditioner BlockDP(offsets);
495 BlockDP.SetDiagonalBlock(0, pc_r.get());
496 BlockDP.SetDiagonalBlock(1, pc_i.get());
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);
520 for (int i = 0; i < Geometry::NumGeom; ++i)
521 {
523 }
524
525 real_t L2Error_Re = x.real().ComputeL2Error(E_ex_Re, irs,
526 pml->GetMarkedPMLElements());
527 real_t L2Error_Im = x.imag().ComputeL2Error(E_ex_Im, irs,
528 pml->GetMarkedPMLElements());
529
530 ComplexGridFunction x_gf0(fespace);
531 x_gf0 = 0.0;
532 real_t 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 real_t t = (real_t)(i % num_frames) / num_frames;
600 ostringstream oss;
601 oss << "Harmonic Solution (t = " << t << " T)";
602
603 add(cos(real_t(2.0 * M_PI) * t), x.real(),
604 sin(real_t(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
620void source(const Vector &x, Vector &f)
621{
622 Vector center(dim);
623 real_t r = 0.0;
624 for (int i = 0; i < dim; ++i)
625 {
626 center(i) = 0.5_r * (comp_domain_bdr(i, 0) + comp_domain_bdr(i, 1));
627 r += pow2(x[i] - center[i]);
628 }
629 real_t n = 5_r * omega * sqrt(epsilon * mu) / real_t(M_PI);
630 real_t coeff = pow2(n) / real_t(M_PI);
631 real_t alpha = -pow2(n) * r;
632 f = 0.0;
633 f[0] = coeff * exp(alpha);
634}
635
636void maxwell_solution(const Vector &x, vector<complex<real_t>> &E)
637{
638 // Initialize
639 for (int i = 0; i < dim; ++i)
640 {
641 E[i] = 0.0;
642 }
643
644 constexpr complex<real_t> zi = complex<real_t>(0., 1.);
645 real_t 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 real_t x0 = x(0) + shift(0);
661 real_t x1 = x(1) + shift(1);
662 real_t r = sqrt(x0 * x0 + x1 * x1);
663 real_t beta = k * r;
664
665 // Bessel functions
666 complex<real_t> Ho, Ho_r, Ho_rr;
667 Ho = real_t(jn(0, beta)) + zi * real_t(yn(0, beta));
668 Ho_r = -k * (real_t(jn(1, beta)) + zi * real_t(yn(1, beta)));
669 Ho_rr = -k * k * (1_r / beta *
670 (real_t(jn(1, beta)) + zi * real_t(yn(1, beta))) -
671 (real_t(jn(2, beta)) + zi * real_t(yn(2, beta))));
672
673 // First derivatives
674 real_t r_x = x0 / r;
675 real_t r_y = x1 / r;
676 real_t r_xy = -(r_x / r) * r_y;
677 real_t r_xx = (1_r / r) * (1_r - r_x * r_x);
678
679 complex<real_t> val, val_xx, val_xy;
680 val = real_t(0.25) * zi * Ho;
681 val_xx = real_t(0.25) * zi * (r_xx * Ho_r + r_x * r_x * Ho_rr);
682 val_xy = real_t(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 real_t x0 = x(0) + shift(0);
689 real_t x1 = x(1) + shift(1);
690 real_t x2 = x(2) + shift(2);
691 real_t r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
692
693 real_t r_x = x0 / r;
694 real_t r_y = x1 / r;
695 real_t r_z = x2 / r;
696 real_t r_xx = (1_r / r) * (1_r - r_x * r_x);
697 real_t r_yx = -(r_y / r) * r_x;
698 real_t r_zx = -(r_z / r) * r_x;
699
700 complex<real_t> val, val_r, val_rr;
701 val = exp(zi * k * r) / r;
702 val_r = val / r * (zi * k * r - 1_r);
703 val_rr = val / (r * r) * (-k * k * r * r
704 - real_t(2) * zi * k * r + real_t(2));
705
706 complex<real_t> 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<real_t> alpha = zi * k / real_t(4) / (real_t) 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 real_t k10 = sqrt(k * k - real_t(M_PI * M_PI));
724 E[1] = -zi * k / (real_t) M_PI *
725 sin((real_t) M_PI*x(2))*exp(zi * k10 * x(0));
726 }
727 else if (dim == 2)
728 {
729 E[1] = -zi * k / (real_t) M_PI * exp(zi * k * x(0));
730 }
731 break;
732 }
733 default:
734 break;
735 }
736}
737
738void E_exact_Re(const Vector &x, Vector &E)
739{
740 vector<complex<real_t>> Eval(E.Size());
741 maxwell_solution(x, Eval);
742 for (int i = 0; i < dim; ++i)
743 {
744 E[i] = Eval[i].real();
745 }
746}
747
748void E_exact_Im(const Vector &x, Vector &E)
749{
750 vector<complex<real_t>> Eval(E.Size());
751 maxwell_solution(x, Eval);
752 for (int i = 0; i < dim; ++i)
753 {
754 E[i] = Eval[i].imag();
755 }
756}
757
758void E_bdr_data_Re(const Vector &x, Vector &E)
759{
760 E = 0.0;
761 bool in_pml = false;
762
763 for (int i = 0; i < dim; ++i)
764 {
765 // check if in PML
766 if (x(i) - comp_domain_bdr(i, 0) < 0_r ||
767 x(i) - comp_domain_bdr(i, 1) > 0_r)
768 {
769 in_pml = true;
770 break;
771 }
772 }
773 if (!in_pml)
774 {
775 vector<complex<real_t>> Eval(E.Size());
776 maxwell_solution(x, Eval);
777 for (int i = 0; i < dim; ++i)
778 {
779 E[i] = Eval[i].real();
780 }
781 }
782}
783
784// Define bdr_data solution
785void E_bdr_data_Im(const Vector &x, Vector &E)
786{
787 E = 0.0;
788 bool in_pml = false;
789
790 for (int i = 0; i < dim; ++i)
791 {
792 // check if in PML
793 if (x(i) - comp_domain_bdr(i, 0) < 0_r ||
794 x(i) - comp_domain_bdr(i, 1) > 0_r)
795 {
796 in_pml = true;
797 break;
798 }
799 }
800 if (!in_pml)
801 {
802 vector<complex<real_t>> Eval(E.Size());
803 maxwell_solution(x, Eval);
804 for (int i = 0; i < dim; ++i)
805 {
806 E[i] = Eval[i].imag();
807 }
808 }
809}
810
811void detJ_JT_J_inv_Re(const Vector &x, PML * pml, Vector &D)
812{
813 vector<complex<real_t>> dxs(dim);
814 complex<real_t> det(1.0, 0.0);
815 pml->StretchFunction(x, dxs);
816
817 for (int i = 0; i < dim; ++i)
818 {
819 det *= dxs[i];
820 }
821
822 for (int i = 0; i < dim; ++i)
823 {
824 D(i) = (det / pow2(dxs[i])).real();
825 }
826}
827
828void detJ_JT_J_inv_Im(const Vector &x, PML * pml, Vector &D)
829{
830 vector<complex<real_t>> dxs(dim);
831 complex<real_t> det = 1.0;
832 pml->StretchFunction(x, dxs);
833
834 for (int i = 0; i < dim; ++i)
835 {
836 det *= dxs[i];
837 }
838
839 for (int i = 0; i < dim; ++i)
840 {
841 D(i) = (det / pow2(dxs[i])).imag();
842 }
843}
844
845void detJ_JT_J_inv_abs(const Vector &x, PML * pml, Vector &D)
846{
847 vector<complex<real_t>> dxs(dim);
848 complex<real_t> det = 1.0;
849 pml->StretchFunction(x, dxs);
850
851 for (int i = 0; i < dim; ++i)
852 {
853 det *= dxs[i];
854 }
855
856 for (int i = 0; i < dim; ++i)
857 {
858 D(i) = abs(det / pow2(dxs[i]));
859 }
860}
861
862void detJ_inv_JT_J_Re(const Vector &x, PML * pml, Vector &D)
863{
864 vector<complex<real_t>> dxs(dim);
865 complex<real_t> det(1.0, 0.0);
866 pml->StretchFunction(x, dxs);
867
868 for (int i = 0; i < dim; ++i)
869 {
870 det *= dxs[i];
871 }
872
873 // in the 2D case the coefficient is scalar 1/det(J)
874 if (dim == 2)
875 {
876 D = (1_r / det).real();
877 }
878 else
879 {
880 for (int i = 0; i < dim; ++i)
881 {
882 D(i) = (pow2(dxs[i]) / det).real();
883 }
884 }
885}
886
887void detJ_inv_JT_J_Im(const Vector &x, PML * pml, Vector &D)
888{
889 vector<complex<real_t>> dxs(dim);
890 complex<real_t> det = 1.0;
891 pml->StretchFunction(x, dxs);
892
893 for (int i = 0; i < dim; ++i)
894 {
895 det *= dxs[i];
896 }
897
898 if (dim == 2)
899 {
900 D = (1_r / det).imag();
901 }
902 else
903 {
904 for (int i = 0; i < dim; ++i)
905 {
906 D(i) = (pow2(dxs[i]) / det).imag();
907 }
908 }
909}
910
911void detJ_inv_JT_J_abs(const Vector &x, PML * pml, Vector &D)
912{
913 vector<complex<real_t>> dxs(dim);
914 complex<real_t> det = 1.0;
915 pml->StretchFunction(x, dxs);
916
917 for (int i = 0; i < dim; ++i)
918 {
919 det *= dxs[i];
920 }
921
922 if (dim == 2)
923 {
924 D = abs(1_r / det);
925 }
926 else
927 {
928 for (int i = 0; i < dim; ++i)
929 {
930 D(i) = abs(pow2(dxs[i]) / det);
931 }
932 }
933}
934
935PML::PML(Mesh *mesh_, Array2D<real_t> length_)
936 : mesh(mesh_), length(length_)
937{
938 dim = mesh->Dimension();
939 SetBoundaries();
940}
941
942void PML::SetBoundaries()
943{
944 comp_dom_bdr.SetSize(dim, 2);
945 dom_bdr.SetSize(dim, 2);
946 Vector pmin, pmax;
947 mesh->GetBoundingBox(pmin, pmax);
948 for (int i = 0; i < dim; i++)
949 {
950 dom_bdr(i, 0) = pmin(i);
951 dom_bdr(i, 1) = pmax(i);
952 comp_dom_bdr(i, 0) = dom_bdr(i, 0) + length(i, 0);
953 comp_dom_bdr(i, 1) = dom_bdr(i, 1) - length(i, 1);
954 }
955}
956
957void PML::SetAttributes(Mesh *mesh_)
958{
959 // Initialize bdr attributes
960 for (int i = 0; i < mesh_->GetNBE(); ++i)
961 {
962 mesh_->GetBdrElement(i)->SetAttribute(i+1);
963 }
964
965 int nrelem = mesh_->GetNE();
966
967 elems.SetSize(nrelem);
968
969 // Loop through the elements and identify which of them are in the PML
970 for (int i = 0; i < nrelem; ++i)
971 {
972 elems[i] = 1;
973 bool in_pml = false;
974 Element *el = mesh_->GetElement(i);
975 Array<int> vertices;
976
977 // Initialize attribute
978 el->SetAttribute(1);
979 el->GetVertices(vertices);
980 int nrvert = vertices.Size();
981
982 // Check if any vertex is in the PML
983 for (int iv = 0; iv < nrvert; ++iv)
984 {
985 int vert_idx = vertices[iv];
986 real_t *coords = mesh_->GetVertex(vert_idx);
987 for (int comp = 0; comp < dim; ++comp)
988 {
989 if (coords[comp] > comp_dom_bdr(comp, 1) ||
990 coords[comp] < comp_dom_bdr(comp, 0))
991 {
992 in_pml = true;
993 break;
994 }
995 }
996 }
997 if (in_pml)
998 {
999 elems[i] = 0;
1000 el->SetAttribute(2);
1001 }
1002 }
1003 mesh_->SetAttributes();
1004}
1005
1006void PML::StretchFunction(const Vector &x,
1007 vector<complex<real_t>> &dxs)
1008{
1009 constexpr complex<real_t> zi = complex<real_t>(0., 1.);
1010
1011 real_t n = 2.0;
1012 real_t c = 5.0;
1013 real_t coeff;
1014 real_t k = omega * sqrt(epsilon * mu);
1015
1016 // Stretch in each direction independently
1017 for (int i = 0; i < dim; ++i)
1018 {
1019 dxs[i] = 1.0;
1020 if (x(i) >= comp_domain_bdr(i, 1))
1021 {
1022 coeff = n * c / k / pow(length(i, 1), n);
1023 dxs[i] = 1_r + zi * coeff *
1024 abs(pow(x(i) - comp_domain_bdr(i, 1), n - 1_r));
1025 }
1026 if (x(i) <= comp_domain_bdr(i, 0))
1027 {
1028 coeff = n * c / k / pow(length(i, 0), n);
1029 dxs[i] = 1_r + zi * coeff *
1030 abs(pow(x(i) - comp_domain_bdr(i, 0), n - 1_r));
1031 }
1032 }
1033}
Dynamic 2D array using row-major layout.
Definition array.hpp:372
void SetSize(int m, int n)
Definition array.hpp:385
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets Operator::DiagonalPolicy used upon construction of the linear system. Policies include:
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
A class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff, VectorCoefficient &imag_coeff, Array< int > &attr)
@ HERMITIAN
Native convention for Hermitian operators.
@ BLOCK_SYMMETRIC
Alternate convention for damping operators.
Specialization of the ComplexOperator built from a pair of Sparse Matrices.
Interface with UMFPack solver specialized for ComplexSparseMatrix This approach avoids forming a mono...
virtual void Mult(const Vector &b, Vector &x) const
This is solving the system A x = b.
real_t Control[UMFPACK_CONTROL]
void SetPrintLevel(int print_lvl)
A coefficient that is constant across space and time.
Integrator for for Nedelec elements.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Abstract data type element.
Definition element.hpp:29
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
Definition element.hpp:58
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:588
GMRES method.
Definition solvers.hpp:547
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition solvers.hpp:559
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the GMRES method.
Definition solvers.cpp:980
Data type for Gauss-Seidel smoother of sparse matrix.
static const int NumGeom
Definition geom.hpp:42
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
Definition geom.hpp:71
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:179
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
Geometry::Type GetBdrElementGeometry(int i) const
Definition mesh.hpp:1376
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition mesh.hpp:1339
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
void GetBoundingBox(Vector &min, Vector &max, int ref=2)
Returns the minimum and maximum corners of the mesh bounding box.
Definition mesh.cpp:138
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1298
ElementTransformation * GetBdrElementTransformation(int i)
Returns a pointer to the transformation defining the i-th boundary element.
Definition mesh.cpp:506
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition mesh.cpp:1604
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1265
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition fe_coll.hpp:465
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:313
@ DIAG_ONE
Set the diagonal value to one.
Definition operator.hpp:50
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
Vector coefficient defined as a product of scalar and vector coefficients.
Scaled Operator B: x -> a A(x).
Definition operator.hpp:727
Data type sparse matrix.
Definition sparsemat.hpp:51
Base class for vector Coefficients that optionally depend on time and space.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
A general vector function coefficient.
Derived vector coefficient that has the value of the parent vector where it is active and is zero oth...
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
Vector beta
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
real_t freq
Definition ex24.cpp:54
real_t omega
Definition ex25.cpp:142
void detJ_JT_J_inv_abs(const Vector &x, PML *pml, Vector &D)
Definition ex25.cpp:845
void detJ_inv_JT_J_Im(const Vector &x, PML *pml, Vector &D)
Definition ex25.cpp:887
T pow2(const T &x)
Definition ex25.cpp:146
void E_bdr_data_Re(const Vector &x, Vector &E)
Definition ex25.cpp:758
bool exact_known
Definition ex25.cpp:144
void detJ_JT_J_inv_Re(const Vector &x, PML *pml, Vector &D)
Definition ex25.cpp:811
int dim
Definition ex25.cpp:143
real_t mu
Definition ex25.cpp:140
real_t epsilon
Definition ex25.cpp:141
void detJ_inv_JT_J_abs(const Vector &x, PML *pml, Vector &D)
Definition ex25.cpp:911
void maxwell_solution(const Vector &x, vector< complex< real_t > > &Eval)
Definition ex25.cpp:636
void source(const Vector &x, Vector &f)
Definition ex25.cpp:620
Array2D< real_t > comp_domain_bdr
Definition ex25.cpp:137
void E_exact_Re(const Vector &x, Vector &E)
Definition ex25.cpp:738
void detJ_JT_J_inv_Im(const Vector &x, PML *pml, Vector &D)
Definition ex25.cpp:828
void detJ_inv_JT_J_Re(const Vector &x, PML *pml, Vector &D)
Definition ex25.cpp:862
prob_type prob
Definition ex25.cpp:156
prob_type
Definition ex25.cpp:149
Definition ex25.cpp:154
@ disc
Definition ex25.cpp:151
@ lshape
Definition ex25.cpp:152
@ beam
Definition ex25.cpp:150
@ fichera
Definition ex25.cpp:153
void E_exact_Im(const Vector &x, Vector &E)
Definition ex25.cpp:748
Array2D< real_t > domain_bdr
Definition ex25.cpp:138
void E_bdr_data_Im(const Vector &x, Vector &E)
Definition ex25.cpp:785
real_t omega
Definition ex25p.cpp:141
int dim
Definition ex25p.cpp:142
real_t mu
Definition ex25p.cpp:139
real_t epsilon
Definition ex25p.cpp:140
Array2D< real_t > comp_domain_bdr
Definition ex25p.cpp:136
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
Geometry Geometries
Definition fe.cpp:49
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:316
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30