MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex25p.cpp
Go to the documentation of this file.
1// MFEM Example 25 - Parallel Version
2//
3// Compile with: make ex25p
4//
5// Sample runs: mpirun -np 4 ex25p -o 2 -f 1.0 -rs 1 -rp 1 -prob 0
6// mpirun -np 4 ex25p -o 3 -f 10.0 -rs 1 -rp 1 -prob 1
7// mpirun -np 4 ex25p -o 3 -f 5.0 -rs 3 -rp 1 -prob 2
8// mpirun -np 4 ex25p -o 2 -f 1.0 -rs 1 -rp 1 -prob 3
9// mpirun -np 4 ex25p -o 2 -f 1.0 -rs 2 -rp 2 -prob 0 -m ../data/beam-quad.mesh
10// mpirun -np 4 ex25p -o 2 -f 8.0 -rs 2 -rp 2 -prob 4 -m ../data/inline-quad.mesh
11// mpirun -np 4 ex25p -o 2 -f 2.0 -rs 1 -rp 1 -prob 4 -m ../data/inline-hex.mesh
12//
13// Device sample runs:
14// mpirun -np 4 ex25p -o 1 -f 3.0 -rs 3 -rp 1 -prob 2 -pa -d cuda
15// mpirun -np 4 ex25p -o 2 -f 1.0 -rs 1 -rp 1 -prob 3 -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
43using namespace std;
44using namespace mfem;
45
46// Class for setting up a simple Cartesian PML region
47class PML
48{
49private:
50 Mesh *mesh;
51
52 int dim;
53
54 // Length of the PML Region in each direction
55 Array2D<real_t> length;
56
57 // Computational Domain Boundary
58 Array2D<real_t> comp_dom_bdr;
59
60 // Domain Boundary
61 Array2D<real_t> 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
70public:
71 // Constructor
72 PML(Mesh *mesh_,Array2D<real_t> length_);
73
74 // Return Computational Domain Boundary
75 Array2D<real_t> GetCompDomainBdr() {return comp_dom_bdr;}
76
77 // Return Domain Boundary
78 Array2D<real_t> 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(ParMesh *pmesh);
85
86 // PML complex stretching function
87 void StretchFunction(const Vector &x, vector<complex<real_t>> &dxs);
88};
89
90// Class for returning the PML coefficients of the bilinear form
91class PMLDiagMatrixCoefficient : public VectorCoefficient
92{
93private:
94 PML * pml = nullptr;
95 void (*Function)(const Vector &, PML *, Vector &);
96public:
97 PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, PML *,
98 Vector &),
99 PML * pml_)
100 : VectorCoefficient(dim), pml(pml_), Function(F)
101 {}
102
104
105 virtual void Eval(Vector &K, ElementTransformation &T,
106 const IntegrationPoint &ip)
107 {
108 real_t x[3];
109 Vector transip(x, 3);
110 T.Transform(ip, transip);
111 K.SetSize(vdim);
112 (*Function)(transip, pml, K);
113 }
114};
115
116void maxwell_solution(const Vector &x, vector<complex<real_t>> &Eval);
117
118void E_bdr_data_Re(const Vector &x, Vector &E);
119void E_bdr_data_Im(const Vector &x, Vector &E);
120
121void E_exact_Re(const Vector &x, Vector &E);
122void E_exact_Im(const Vector &x, Vector &E);
123
124void 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
128void detJ_JT_J_inv_Re(const Vector &x, PML * pml, Vector & D);
129void detJ_JT_J_inv_Im(const Vector &x, PML * pml, Vector & D);
130void detJ_JT_J_inv_abs(const Vector &x, PML * pml, Vector & D);
131
132void detJ_inv_JT_J_Re(const Vector &x, PML * pml, Vector & D);
133void detJ_inv_JT_J_Im(const Vector &x, PML * pml, Vector & D);
134void detJ_inv_JT_J_abs(const Vector &x, PML * pml, Vector & D);
135
138
139real_t mu = 1.0;
142int dim;
143bool exact_known = false;
144
145template <typename T> T pow2(const T &x) { return x*x; }
146
148{
149 beam, // Wave propagating in a beam-like domain
150 disc, // Point source propagating in the square-disc domain
151 lshape, // Point source propagating in the L-shape domain
152 fichera, // Point source propagating in the fichera domain
153 load_src // Approximated point source with PML all around
156
157int main(int argc, char *argv[])
158{
159 // 1. Initialize MPI and HYPRE.
160 Mpi::Init(argc, argv);
161 int num_procs = Mpi::WorldSize();
162 int myid = Mpi::WorldRank();
163 Hypre::Init();
164
165 // 2. Parse command-line options.
166 const char *mesh_file = nullptr;
167 int order = 1;
168 int ref_levels = 1;
169 int par_ref_levels = 2;
170 int iprob = 4;
171 real_t freq = 5.0;
172 bool herm_conv = true;
173 bool slu_solver = false;
174 bool mumps_solver = false;
175 bool strumpack_solver = false;
176 bool visualization = 1;
177 bool pa = false;
178 const char *device_config = "cpu";
179
180 OptionsParser args(argc, argv);
181 args.AddOption(&mesh_file, "-m", "--mesh",
182 "Mesh file to use.");
183 args.AddOption(&order, "-o", "--order",
184 "Finite element order (polynomial degree).");
185 args.AddOption(&iprob, "-prob", "--problem", "Problem case"
186 " 0: beam, 1: disc, 2: lshape, 3: fichera, 4: General");
187 args.AddOption(&ref_levels, "-rs", "--refinements-serial",
188 "Number of serial refinements");
189 args.AddOption(&par_ref_levels, "-rp", "--refinements-parallel",
190 "Number of parallel refinements");
191 args.AddOption(&mu, "-mu", "--permeability",
192 "Permeability of free space (or 1/(spring constant)).");
193 args.AddOption(&epsilon, "-eps", "--permittivity",
194 "Permittivity of free space (or mass constant).");
195 args.AddOption(&freq, "-f", "--frequency",
196 "Frequency (in Hz).");
197 args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
198 "--no-hermitian", "Use convention for Hermitian operators.");
199#ifdef MFEM_USE_SUPERLU
200 args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
201 "--no-superlu", "Use the SuperLU Solver.");
202#endif
203#ifdef MFEM_USE_MUMPS
204 args.AddOption(&mumps_solver, "-mumps", "--mumps-solver", "-no-mumps",
205 "--no-mumps-solver", "Use the MUMPS Solver.");
206#endif
207#ifdef MFEM_USE_STRUMPACK
208 args.AddOption(&strumpack_solver, "-strumpack", "--strumpack-solver",
209 "-no-strumpack", "--no-strumpack-solver",
210 "Use the STRUMPACK Solver.");
211#endif
212 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
213 "--no-visualization",
214 "Enable or disable GLVis visualization.");
215 args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
216 "--no-partial-assembly", "Enable Partial Assembly.");
217 args.AddOption(&device_config, "-d", "--device",
218 "Device configuration string, see Device::Configure().");
219 args.Parse();
220 if (slu_solver + mumps_solver + strumpack_solver > 1)
221 {
222 if (myid == 0)
223 cout << "WARNING: More than one of SuperLU, MUMPS, and STRUMPACK have"
224 << " been selected, please choose only one." << endl
225 << " Defaulting to SuperLU." << endl;
226 mumps_solver = false;
227 strumpack_solver = false;
228 }
229
230 if (iprob > 4) { iprob = 4; }
231 prob = (prob_type)iprob;
232
233 // 3. Enable hardware devices such as GPUs, and programming models such as
234 // CUDA, OCCA, RAJA and OpenMP based on command line options.
235 Device device(device_config);
236 if (myid == 0) { device.Print(); }
237
238 // 4. Setup the (serial) mesh on all processors.
239 if (!mesh_file)
240 {
241 exact_known = true;
242 switch (prob)
243 {
244 case beam:
245 mesh_file = "../data/beam-hex.mesh";
246 break;
247 case disc:
248 mesh_file = "../data/square-disc.mesh";
249 break;
250 case lshape:
251 mesh_file = "../data/l-shape.mesh";
252 break;
253 case fichera:
254 mesh_file = "../data/fichera.mesh";
255 break;
256 default:
257 exact_known = false;
258 mesh_file = "../data/inline-quad.mesh";
259 break;
260 }
261 }
262
263 if (!args.Good())
264 {
265 if (myid == 0)
266 {
267 args.PrintUsage(cout);
268 }
269 return 1;
270 }
271 if (myid == 0)
272 {
273 args.PrintOptions(cout);
274 }
275
276 Mesh * mesh = new Mesh(mesh_file, 1, 1);
277 dim = mesh->Dimension();
278
279 // Angular frequency
280 omega = real_t(2.0 * M_PI) * freq;
281
282 // Setup PML length
283 Array2D<real_t> length(dim, 2); length = 0.0;
284
285 // 5. Setup the Cartesian PML region.
286 switch (prob)
287 {
288 case disc:
289 length = 0.2;
290 break;
291 case lshape:
292 length(0, 0) = 0.1;
293 length(1, 0) = 0.1;
294 break;
295 case fichera:
296 length(0, 1) = 0.5;
297 length(1, 1) = 0.5;
298 length(2, 1) = 0.5;
299 break;
300 case beam:
301 length(0, 1) = 2.0;
302 break;
303 default:
304 length = 0.25;
305 break;
306 }
307 PML * pml = new PML(mesh,length);
308 comp_domain_bdr = pml->GetCompDomainBdr();
309 domain_bdr = pml->GetDomainBdr();
310
311 // 6. Refine the serial mesh on all processors to increase the resolution.
312 for (int l = 0; l < ref_levels; l++)
313 {
314 mesh->UniformRefinement();
315 }
316
317 // 7. Define a parallel mesh by a partitioning of the serial mesh.
318 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
319 delete mesh;
320 {
321 for (int l = 0; l < par_ref_levels; l++)
322 {
323 pmesh->UniformRefinement();
324 }
325 }
326
327 // 8. Set element attributes in order to distinguish elements in the PML
328 pml->SetAttributes(pmesh);
329
330 // 9. Define a parallel finite element space on the parallel mesh. Here we
331 // use the Nedelec finite elements of the specified order.
333 ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
334 HYPRE_BigInt size = fespace->GlobalTrueVSize();
335 if (myid == 0)
336 {
337 cout << "Number of finite element unknowns: " << size << endl;
338 }
339
340 // 10. Determine the list of true (i.e. parallel conforming) essential
341 // boundary dofs. In this example, the boundary conditions are defined
342 // based on the specific mesh and the problem type.
343 Array<int> ess_tdof_list;
344 Array<int> ess_bdr;
345 if (pmesh->bdr_attributes.Size())
346 {
347 ess_bdr.SetSize(pmesh->bdr_attributes.Max());
348 ess_bdr = 1;
349 if (prob == lshape || prob == fichera)
350 {
351 ess_bdr = 0;
352 for (int j = 0; j < pmesh->GetNBE(); j++)
353 {
354 Vector center(dim);
355 int bdrgeom = pmesh->GetBdrElementGeometry(j);
357 tr->Transform(Geometries.GetCenter(bdrgeom),center);
358 int k = pmesh->GetBdrAttribute(j);
359 switch (prob)
360 {
361 case lshape:
362 if (center[0] == 1_r || center[0] == 0.5_r ||
363 center[1] == 0.5_r)
364 {
365 ess_bdr[k - 1] = 1;
366 }
367 break;
368 case fichera:
369 if (center[0] == -1_r || center[0] == 0_r ||
370 center[1] == 0_r || center[2] == 0_r)
371 {
372 ess_bdr[k - 1] = 1;
373 }
374 break;
375 default:
376 break;
377 }
378 }
379 }
380 }
381 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
382
383 // 11. Setup Complex Operator convention
386
387 // 12. Set up the parallel linear form b(.) which corresponds to the
388 // right-hand side of the FEM linear system.
390 ParComplexLinearForm b(fespace, conv);
391 if (prob == load_src)
392 {
393 b.AddDomainIntegrator(NULL, new VectorFEDomainLFIntegrator(f));
394 }
395 b.Vector::operator=(0.0);
396 b.Assemble();
397
398 // 13. Define the solution vector x as a parallel complex finite element grid
399 // function corresponding to fespace.
400 ParComplexGridFunction x(fespace);
401 x = 0.0;
404 x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
405
406 // 14. Set up the parallel sesquilinear form a(.,.)
407 //
408 // In Comp
409 // Domain: 1/mu (Curl E, Curl F) - omega^2 * epsilon (E,F)
410 //
411 // In PML: 1/mu (1/det(J) J^T J Curl E, Curl F)
412 // - omega^2 * epsilon (det(J) * (J^T J)^-1 * E, F)
413 //
414 // where J denotes the Jacobian Matrix of the PML Stretching function
415 Array<int> attr;
416 Array<int> attrPML;
417 if (pmesh->attributes.Size())
418 {
419 attr.SetSize(pmesh->attributes.Max());
420 attrPML.SetSize(pmesh->attributes.Max());
421 attr = 0; attr[0] = 1;
422 attrPML = 0;
423 if (pmesh->attributes.Max() > 1)
424 {
425 attrPML[1] = 1;
426 }
427 }
428
429 ConstantCoefficient muinv(1_r / mu);
431 RestrictedCoefficient restr_muinv(muinv,attr);
432 RestrictedCoefficient restr_omeg(omeg,attr);
433
434 // Integrators inside the computational domain (excluding the PML region)
435 ParSesquilinearForm a(fespace, conv);
438
439 int cdim = (dim == 2) ? 1 : dim;
440 PMLDiagMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, pml);
441 PMLDiagMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, pml);
442 ScalarVectorProductCoefficient c1_Re(muinv,pml_c1_Re);
443 ScalarVectorProductCoefficient c1_Im(muinv,pml_c1_Im);
444 VectorRestrictedCoefficient restr_c1_Re(c1_Re,attrPML);
445 VectorRestrictedCoefficient restr_c1_Im(c1_Im,attrPML);
446
447 PMLDiagMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,pml);
448 PMLDiagMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,pml);
449 ScalarVectorProductCoefficient c2_Re(omeg,pml_c2_Re);
450 ScalarVectorProductCoefficient c2_Im(omeg,pml_c2_Im);
451 VectorRestrictedCoefficient restr_c2_Re(c2_Re,attrPML);
452 VectorRestrictedCoefficient restr_c2_Im(c2_Im,attrPML);
453
454 // Integrators inside the PML region
456 new CurlCurlIntegrator(restr_c1_Im));
458 new VectorFEMassIntegrator(restr_c2_Im));
459
460 // 15. Assemble the parallel bilinear form and the corresponding linear
461 // system, applying any necessary transformations such as: parallel
462 // assembly, eliminating boundary conditions, applying conforming
463 // constraints for non-conforming AMR, etc.
464 if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
465 a.Assemble();
466
467 OperatorPtr Ah;
468 Vector B, X;
469 a.FormLinearSystem(ess_tdof_list, x, b, Ah, X, B);
470
471 // 16. Solve using a direct or an iterative solver
472#ifdef MFEM_USE_SUPERLU
473 if (!pa && slu_solver)
474 {
475 // Transform to monolithic HypreParMatrix
476 HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
477 SuperLURowLocMatrix SA(*A);
478 SuperLUSolver superlu(MPI_COMM_WORLD);
479 superlu.SetPrintStatistics(false);
480 superlu.SetSymmetricPattern(false);
482 superlu.SetOperator(SA);
483 superlu.Mult(B, X);
484 delete A;
485 }
486#endif
487#ifdef MFEM_USE_STRUMPACK
488 if (!pa && strumpack_solver)
489 {
490 HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
492 STRUMPACKSolver strumpack(MPI_COMM_WORLD, argc, argv);
493 strumpack.SetPrintFactorStatistics(false);
494 strumpack.SetPrintSolveStatistics(false);
495 strumpack.SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
496 strumpack.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
497 strumpack.SetMatching(strumpack::MatchingJob::NONE);
498 strumpack.SetCompression(strumpack::CompressionType::NONE);
499 strumpack.SetFromCommandLine();
500 strumpack.SetOperator(SA);
501 strumpack.Mult(B, X);
502 delete A;
503 }
504#endif
505#ifdef MFEM_USE_MUMPS
506 if (!pa && mumps_solver)
507 {
508 HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
509 MUMPSSolver mumps(A->GetComm());
510 mumps.SetPrintLevel(0);
511 mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
512 mumps.SetOperator(*A);
513 mumps.Mult(B, X);
514 delete A;
515 }
516#endif
517 // 16a. Set up the parallel Bilinear form a(.,.) for the preconditioner
518 //
519 // In Comp
520 // Domain: 1/mu (Curl E, Curl F) + omega^2 * epsilon (E,F)
521 //
522 // In PML: 1/mu (abs(1/det(J) J^T J) Curl E, Curl F)
523 // + omega^2 * epsilon (abs(det(J) * (J^T J)^-1) * E, F)
524 if (pa || (!slu_solver && !mumps_solver && !strumpack_solver))
525 {
527 RestrictedCoefficient restr_absomeg(absomeg,attr);
528
529 ParBilinearForm prec(fespace);
532
533 PMLDiagMatrixCoefficient pml_c1_abs(cdim,detJ_inv_JT_J_abs, pml);
534 ScalarVectorProductCoefficient c1_abs(muinv,pml_c1_abs);
535 VectorRestrictedCoefficient restr_c1_abs(c1_abs,attrPML);
536
537 PMLDiagMatrixCoefficient pml_c2_abs(dim, detJ_JT_J_inv_abs,pml);
538 ScalarVectorProductCoefficient c2_abs(absomeg,pml_c2_abs);
539 VectorRestrictedCoefficient restr_c2_abs(c2_abs,attrPML);
540
543
544 if (pa) { prec.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
545 prec.Assemble();
546
547 // 16b. Define and apply a parallel GMRES solver for AU=B with a block
548 // diagonal preconditioner based on hypre's AMS preconditioner.
549 Array<int> offsets(3);
550 offsets[0] = 0;
551 offsets[1] = fespace->GetTrueVSize();
552 offsets[2] = fespace->GetTrueVSize();
553 offsets.PartialSum();
554
555 std::unique_ptr<Operator> pc_r;
556 std::unique_ptr<Operator> pc_i;
557 int s = (conv == ComplexOperator::HERMITIAN) ? -1 : 1;
558 if (pa)
559 {
560 // Jacobi Smoother
561 pc_r.reset(new OperatorJacobiSmoother(prec, ess_tdof_list));
562 pc_i.reset(new ScaledOperator(pc_r.get(), s));
563 }
564 else
565 {
566 OperatorPtr PCOpAh;
567 prec.FormSystemMatrix(ess_tdof_list, PCOpAh);
568
569 // Hypre AMS
570 pc_r.reset(new HypreAMS(*PCOpAh.As<HypreParMatrix>(), fespace));
571 pc_i.reset(new ScaledOperator(pc_r.get(), s));
572 }
573
574 BlockDiagonalPreconditioner BlockDP(offsets);
575 BlockDP.SetDiagonalBlock(0, pc_r.get());
576 BlockDP.SetDiagonalBlock(1, pc_i.get());
577
578 GMRESSolver gmres(MPI_COMM_WORLD);
579 gmres.SetPrintLevel(1);
580 gmres.SetKDim(200);
581 gmres.SetMaxIter(pa ? 5000 : 2000);
582 gmres.SetRelTol(1e-5);
583 gmres.SetAbsTol(0.0);
584 gmres.SetOperator(*Ah);
585 gmres.SetPreconditioner(BlockDP);
586 gmres.Mult(B, X);
587 }
588
589 // 17. Recover the parallel grid function corresponding to X. This is the
590 // local finite element solution on each processor.
591 a.RecoverFEMSolution(X, b, x);
592
593 // If exact is known compute the error
594 if (exact_known)
595 {
598 int order_quad = max(2, 2 * order + 1);
600 for (int i = 0; i < Geometry::NumGeom; ++i)
601 {
602 irs[i] = &(IntRules.Get(i, order_quad));
603 }
604
605 real_t L2Error_Re = x.real().ComputeL2Error(E_ex_Re, irs,
606 pml->GetMarkedPMLElements());
607 real_t L2Error_Im = x.imag().ComputeL2Error(E_ex_Im, irs,
608 pml->GetMarkedPMLElements());
609
610 ParComplexGridFunction x_gf0(fespace);
611 x_gf0 = 0.0;
612 real_t norm_E_Re, norm_E_Im;
613 norm_E_Re = x_gf0.real().ComputeL2Error(E_ex_Re, irs,
614 pml->GetMarkedPMLElements());
615 norm_E_Im = x_gf0.imag().ComputeL2Error(E_ex_Im, irs,
616 pml->GetMarkedPMLElements());
617
618 if (myid == 0)
619 {
620 cout << "\n Relative Error (Re part): || E_h - E || / ||E|| = "
621 << L2Error_Re / norm_E_Re
622 << "\n Relative Error (Im part): || E_h - E || / ||E|| = "
623 << L2Error_Im / norm_E_Im
624 << "\n Total Error: "
625 << sqrt(L2Error_Re*L2Error_Re + L2Error_Im*L2Error_Im) << "\n\n";
626 }
627 }
628
629 // 18. Save the refined mesh and the solution in parallel. This output can be
630 // viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
631 {
632 ostringstream mesh_name, sol_r_name, sol_i_name;
633 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
634 sol_r_name << "ex25p-sol_r." << setfill('0') << setw(6) << myid;
635 sol_i_name << "ex25p-sol_i." << setfill('0') << setw(6) << myid;
636
637 ofstream mesh_ofs(mesh_name.str().c_str());
638 mesh_ofs.precision(8);
639 pmesh->Print(mesh_ofs);
640
641 ofstream sol_r_ofs(sol_r_name.str().c_str());
642 ofstream sol_i_ofs(sol_i_name.str().c_str());
643 sol_r_ofs.precision(8);
644 sol_i_ofs.precision(8);
645 x.real().Save(sol_r_ofs);
646 x.imag().Save(sol_i_ofs);
647 }
648
649 // 19. Send the solution by socket to a GLVis server.
650 if (visualization)
651 {
652 // Define visualization keys for GLVis (see GLVis documentation)
653 string keys;
654 keys = (dim == 3) ? "keys macF\n" : keys = "keys amrRljcUUuu\n";
655 if (prob == beam && dim == 3) {keys = "keys macFFiYYYYYYYYYYYYYYYYYY\n";}
656 if (prob == beam && dim == 2) {keys = "keys amrRljcUUuuu\n"; }
657
658 char vishost[] = "localhost";
659 int visport = 19916;
660
661 {
662 socketstream sol_sock_re(vishost, visport);
663 sol_sock_re.precision(8);
664 sol_sock_re << "parallel " << num_procs << " " << myid << "\n"
665 << "solution\n" << *pmesh << x.real() << keys
666 << "window_title 'Solution real part'" << flush;
667 MPI_Barrier(MPI_COMM_WORLD); // try to prevent streams from mixing
668 }
669
670 {
671 socketstream sol_sock_im(vishost, visport);
672 sol_sock_im.precision(8);
673 sol_sock_im << "parallel " << num_procs << " " << myid << "\n"
674 << "solution\n" << *pmesh << x.imag() << keys
675 << "window_title 'Solution imag part'" << flush;
676 MPI_Barrier(MPI_COMM_WORLD); // try to prevent streams from mixing
677 }
678
679 {
680 ParGridFunction x_t(fespace);
681 x_t = x.real();
682
683 socketstream sol_sock(vishost, visport);
684 sol_sock.precision(8);
685 sol_sock << "parallel " << num_procs << " " << myid << "\n"
686 << "solution\n" << *pmesh << x_t << keys << "autoscale off\n"
687 << "window_title 'Harmonic Solution (t = 0.0 T)'"
688 << "pause\n" << flush;
689
690 if (myid == 0)
691 {
692 cout << "GLVis visualization paused."
693 << " Press space (in the GLVis window) to resume it.\n";
694 }
695
696 int num_frames = 32;
697 int i = 0;
698 while (sol_sock)
699 {
700 real_t t = (real_t)(i % num_frames) / num_frames;
701 ostringstream oss;
702 oss << "Harmonic Solution (t = " << t << " T)";
703
705 sin(real_t(2.0*M_PI)*t), x.imag(), x_t);
706 sol_sock << "parallel " << num_procs << " " << myid << "\n";
707 sol_sock << "solution\n" << *pmesh << x_t
708 << "window_title '" << oss.str() << "'" << flush;
709 i++;
710 }
711 }
712 }
713
714 // 20. Free the used memory.
715 delete pml;
716 delete fespace;
717 delete fec;
718 delete pmesh;
719 return 0;
720}
721
722void source(const Vector &x, Vector &f)
723{
724 Vector center(dim);
725 real_t r = 0.0;
726 for (int i = 0; i < dim; ++i)
727 {
728 center(i) = real_t(0.5) * (comp_domain_bdr(i, 0) + comp_domain_bdr(i, 1));
729 r += pow2(x[i] - center[i]);
730 }
731 real_t n = real_t(5) * omega * sqrt(epsilon * mu) / real_t(M_PI);
732 real_t coeff = pow2(n) / real_t(M_PI);
733 real_t alpha = -pow2(n) * r;
734 f = 0.0;
735 f[0] = coeff * exp(alpha);
736}
737
738void maxwell_solution(const Vector &x, vector<complex<real_t>> &E)
739{
740 // Initialize
741 for (int i = 0; i < dim; ++i)
742 {
743 E[i] = 0.0;
744 }
745
746 constexpr complex<real_t> zi = complex<real_t>(0., 1.);
747 real_t k = omega * sqrt(epsilon * mu);
748 switch (prob)
749 {
750 case disc:
751 case lshape:
752 case fichera:
753 {
754 Vector shift(dim);
755 shift = 0.0;
756 if (prob == fichera) { shift = 1.0; }
757 if (prob == disc) { shift = -0.5; }
758 if (prob == lshape) { shift = -1.0; }
759
760 if (dim == 2)
761 {
762 real_t x0 = x(0) + shift(0);
763 real_t x1 = x(1) + shift(1);
764 real_t r = sqrt(x0 * x0 + x1 * x1);
765 real_t beta = k * r;
766
767 // Bessel functions
768 complex<real_t> Ho, Ho_r, Ho_rr;
769 Ho = real_t(jn(0, beta)) + zi * real_t(yn(0, beta));
770 Ho_r = -k * (real_t(jn(1, beta)) + zi * real_t(yn(1, beta)));
771 Ho_rr = -k * k * (1_r / beta *
772 (real_t(jn(1, beta)) + zi * real_t(yn(1, beta))) -
773 (real_t(jn(2, beta)) + zi * real_t(yn(2, beta))));
774
775 // First derivatives
776 real_t r_x = x0 / r;
777 real_t r_y = x1 / r;
778 real_t r_xy = -(r_x / r) * r_y;
779 real_t r_xx = (1_r / r) * (1_r - r_x * r_x);
780
781 complex<real_t> val, val_xx, val_xy;
782 val = real_t(0.25) * zi * Ho;
783 val_xx = real_t(0.25) * zi * (r_xx * Ho_r + r_x * r_x * Ho_rr);
784 val_xy = real_t(0.25) * zi * (r_xy * Ho_r + r_x * r_y * Ho_rr);
785 E[0] = zi / k * (k * k * val + val_xx);
786 E[1] = zi / k * val_xy;
787 }
788 else if (dim == 3)
789 {
790 real_t x0 = x(0) + shift(0);
791 real_t x1 = x(1) + shift(1);
792 real_t x2 = x(2) + shift(2);
793 real_t r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
794
795 real_t r_x = x0 / r;
796 real_t r_y = x1 / r;
797 real_t r_z = x2 / r;
798 real_t r_xx = (1_r / r) * (1_r - r_x * r_x);
799 real_t r_yx = -(r_y / r) * r_x;
800 real_t r_zx = -(r_z / r) * r_x;
801
802 complex<real_t> val, val_r, val_rr;
803 val = exp(zi * k * r) / r;
804 val_r = val / r * (zi * k * r - 1_r);
805 val_rr = val / (r * r) * (-k * k * r * r
806 - real_t(2) * zi * k * r + real_t(2));
807
808 complex<real_t> val_xx, val_yx, val_zx;
809 val_xx = val_rr * r_x * r_x + val_r * r_xx;
810 val_yx = val_rr * r_x * r_y + val_r * r_yx;
811 val_zx = val_rr * r_x * r_z + val_r * r_zx;
812
813 complex<real_t> alpha = zi * k / real_t(4) / (real_t) M_PI / k / k;
814 E[0] = alpha * (k * k * val + val_xx);
815 E[1] = alpha * val_yx;
816 E[2] = alpha * val_zx;
817 }
818 break;
819 }
820 case beam:
821 {
822 // T_10 mode
823 if (dim == 3)
824 {
825 real_t k10 = sqrt(k * k - real_t(M_PI * M_PI));
826 E[1] = -zi * k / (real_t) M_PI *
827 sin((real_t) M_PI*x(2))*exp(zi * k10 * x(0));
828 }
829 else if (dim == 2)
830 {
831 E[1] = -zi * k / (real_t) M_PI * exp(zi * k * x(0));
832 }
833 break;
834 }
835 default:
836 break;
837 }
838}
839
840void E_exact_Re(const Vector &x, Vector &E)
841{
842 vector<complex<real_t>> Eval(E.Size());
843 maxwell_solution(x, Eval);
844 for (int i = 0; i < dim; ++i)
845 {
846 E[i] = Eval[i].real();
847 }
848}
849
850void E_exact_Im(const Vector &x, Vector &E)
851{
852 vector<complex<real_t>> Eval(E.Size());
853 maxwell_solution(x, Eval);
854 for (int i = 0; i < dim; ++i)
855 {
856 E[i] = Eval[i].imag();
857 }
858}
859
860void E_bdr_data_Re(const Vector &x, Vector &E)
861{
862 E = 0.0;
863 bool in_pml = false;
864
865 for (int i = 0; i < dim; ++i)
866 {
867 // check if in PML
868 if (x(i) - comp_domain_bdr(i, 0) < 0_r ||
869 x(i) - comp_domain_bdr(i, 1) > 0_r)
870 {
871 in_pml = true;
872 break;
873 }
874 }
875 if (!in_pml)
876 {
877 vector<complex<real_t>> Eval(E.Size());
878 maxwell_solution(x, Eval);
879 for (int i = 0; i < dim; ++i)
880 {
881 E[i] = Eval[i].real();
882 }
883 }
884}
885
886// Define bdr_data solution
887void E_bdr_data_Im(const Vector &x, Vector &E)
888{
889 E = 0.0;
890 bool in_pml = false;
891
892 for (int i = 0; i < dim; ++i)
893 {
894 // check if in PML
895 if (x(i) - comp_domain_bdr(i, 0) < 0_r ||
896 x(i) - comp_domain_bdr(i, 1) > 0_r)
897 {
898 in_pml = true;
899 break;
900 }
901 }
902 if (!in_pml)
903 {
904 vector<complex<real_t>> Eval(E.Size());
905 maxwell_solution(x, Eval);
906 for (int i = 0; i < dim; ++i)
907 {
908 E[i] = Eval[i].imag();
909 }
910 }
911}
912
913void detJ_JT_J_inv_Re(const Vector &x, PML * pml, Vector & D)
914{
915 vector<complex<real_t>> dxs(dim);
916 complex<real_t> det(1.0, 0.0);
917 pml->StretchFunction(x, dxs);
918
919 for (int i = 0; i < dim; ++i)
920 {
921 det *= dxs[i];
922 }
923
924 for (int i = 0; i < dim; ++i)
925 {
926 D(i) = (det / pow2(dxs[i])).real();
927 }
928}
929
930void detJ_JT_J_inv_Im(const Vector &x, PML * pml, Vector & D)
931{
932 vector<complex<real_t>> dxs(dim);
933 complex<real_t> det = 1.0;
934 pml->StretchFunction(x, dxs);
935
936 for (int i = 0; i < dim; ++i)
937 {
938 det *= dxs[i];
939 }
940
941 for (int i = 0; i < dim; ++i)
942 {
943 D(i) = (det / pow2(dxs[i])).imag();
944 }
945}
946
947void detJ_JT_J_inv_abs(const Vector &x, PML * pml, Vector & D)
948{
949 vector<complex<real_t>> dxs(dim);
950 complex<real_t> det = 1.0;
951 pml->StretchFunction(x, dxs);
952
953 for (int i = 0; i < dim; ++i)
954 {
955 det *= dxs[i];
956 }
957
958 for (int i = 0; i < dim; ++i)
959 {
960 D(i) = abs(det / pow2(dxs[i]));
961 }
962}
963
964void detJ_inv_JT_J_Re(const Vector &x, PML * pml, Vector & D)
965{
966 vector<complex<real_t>> dxs(dim);
967 complex<real_t> det(1.0, 0.0);
968 pml->StretchFunction(x, dxs);
969
970 for (int i = 0; i < dim; ++i)
971 {
972 det *= dxs[i];
973 }
974
975 // in the 2D case the coefficient is scalar 1/det(J)
976 if (dim == 2)
977 {
978 D = (1_r / det).real();
979 }
980 else
981 {
982 for (int i = 0; i < dim; ++i)
983 {
984 D(i) = (pow2(dxs[i]) / det).real();
985 }
986 }
987}
988
989void detJ_inv_JT_J_Im(const Vector &x, PML * pml, Vector & D)
990{
991 vector<complex<real_t>> dxs(dim);
992 complex<real_t> det = 1.0;
993 pml->StretchFunction(x, dxs);
994
995 for (int i = 0; i < dim; ++i)
996 {
997 det *= dxs[i];
998 }
999
1000 if (dim == 2)
1001 {
1002 D = (1_r / det).imag();
1003 }
1004 else
1005 {
1006 for (int i = 0; i < dim; ++i)
1007 {
1008 D(i) = (pow2(dxs[i]) / det).imag();
1009 }
1010 }
1011}
1012
1013void detJ_inv_JT_J_abs(const Vector &x, PML * pml, Vector & D)
1014{
1015 vector<complex<real_t>> dxs(dim);
1016 complex<real_t> det = 1.0;
1017 pml->StretchFunction(x, dxs);
1018
1019 for (int i = 0; i < dim; ++i)
1020 {
1021 det *= dxs[i];
1022 }
1023
1024 if (dim == 2)
1025 {
1026 D = abs(1_r / det);
1027 }
1028 else
1029 {
1030 for (int i = 0; i < dim; ++i)
1031 {
1032 D(i) = abs(pow2(dxs[i]) / det);
1033 }
1034 }
1035}
1036
1037PML::PML(Mesh *mesh_, Array2D<real_t> length_)
1038 : mesh(mesh_), length(length_)
1039{
1040 dim = mesh->Dimension();
1041 SetBoundaries();
1042}
1043
1044void PML::SetBoundaries()
1045{
1046 comp_dom_bdr.SetSize(dim, 2);
1047 dom_bdr.SetSize(dim, 2);
1048 Vector pmin, pmax;
1049 mesh->GetBoundingBox(pmin, pmax);
1050 for (int i = 0; i < dim; i++)
1051 {
1052 dom_bdr(i, 0) = pmin(i);
1053 dom_bdr(i, 1) = pmax(i);
1054 comp_dom_bdr(i, 0) = dom_bdr(i, 0) + length(i, 0);
1055 comp_dom_bdr(i, 1) = dom_bdr(i, 1) - length(i, 1);
1056 }
1057}
1058
1059void PML::SetAttributes(ParMesh *pmesh)
1060{
1061 // Initialize bdr attributes
1062 for (int i = 0; i < pmesh->GetNBE(); ++i)
1063 {
1064 pmesh->GetBdrElement(i)->SetAttribute(i+1);
1065 }
1066
1067 int nrelem = pmesh->GetNE();
1068
1069 // Initialize list with 1
1070 elems.SetSize(nrelem);
1071
1072 // Loop through the elements and identify which of them are in the PML
1073 for (int i = 0; i < nrelem; ++i)
1074 {
1075 elems[i] = 1;
1076 bool in_pml = false;
1077 Element *el = pmesh->GetElement(i);
1078 Array<int> vertices;
1079
1080 // Initialize attribute
1081 el->SetAttribute(1);
1082 el->GetVertices(vertices);
1083 int nrvert = vertices.Size();
1084
1085 // Check if any vertex is in the PML
1086 for (int iv = 0; iv < nrvert; ++iv)
1087 {
1088 int vert_idx = vertices[iv];
1089 real_t *coords = pmesh->GetVertex(vert_idx);
1090 for (int comp = 0; comp < dim; ++comp)
1091 {
1092 if (coords[comp] > comp_dom_bdr(comp, 1) ||
1093 coords[comp] < comp_dom_bdr(comp, 0))
1094 {
1095 in_pml = true;
1096 break;
1097 }
1098 }
1099 }
1100 if (in_pml)
1101 {
1102 elems[i] = 0;
1103 el->SetAttribute(2);
1104 }
1105 }
1106 pmesh->SetAttributes();
1107}
1108
1109void PML::StretchFunction(const Vector &x,
1110 vector<complex<real_t>> &dxs)
1111{
1112 constexpr complex<real_t> zi = complex<real_t>(0., 1.);
1113
1114 real_t n = 2.0;
1115 real_t c = 5.0;
1116 real_t coeff;
1117 real_t k = omega * sqrt(epsilon * mu);
1118
1119 // Stretch in each direction independently
1120 for (int i = 0; i < dim; ++i)
1121 {
1122 dxs[i] = 1.0;
1123 if (x(i) >= comp_domain_bdr(i, 1))
1124 {
1125 coeff = n * c / k / pow(length(i, 1), n);
1126 dxs[i] = 1_r + zi * coeff *
1127 abs(pow(x(i) - comp_domain_bdr(i, 1), n - 1_r));
1128 }
1129 if (x(i) <= comp_domain_bdr(i, 0))
1130 {
1131 coeff = n * c / k / pow(length(i, 0), n);
1132 dxs[i] = 1_r + zi * coeff *
1133 abs(pow(x(i) - comp_domain_bdr(i, 0), n - 1_r));
1134 }
1135 }
1136}
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
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
Adds new Domain Integrator. Assumes ownership of bfi.
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).
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
@ HERMITIAN
Native convention for Hermitian operators.
@ BLOCK_SYMMETRIC
Alternate convention for damping operators.
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
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
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
The Auxiliary-space Maxwell Solver in hypre.
Definition hypre.hpp:1845
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
MPI_Comm GetComm() const
MPI communicator.
Definition hypre.hpp:578
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
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
MUMPS: A Parallel Sparse Direct Solver.
Definition mumps.hpp:39
void SetPrintLevel(int print_lvl)
Set the error print level for MUMPS.
Definition mumps.cpp:445
void Mult(const Vector &x, Vector &y) const
Solve .
Definition mumps.cpp:353
void SetMatrixSymType(MatType mtype)
Set the matrix type.
Definition mumps.cpp:450
void SetOperator(const Operator &op)
Set the Operator and perform factorization.
Definition mumps.cpp:103
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
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
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1265
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
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
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.
Class for parallel bilinear form.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff, VectorCoefficient &imag_coeff, Array< int > &attr)
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
Class for parallel meshes.
Definition pmesh.hpp:34
void SetAttributes() override
Determine the sets of unique attribute values in domain and boundary elements.
Definition pmesh.cpp:1589
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
void SetOperator(const Operator &op)
Set the operator/matrix.
void SetMatching(strumpack::MatchingJob job)
Configure static pivoting for stability.
void SetCompression(strumpack::CompressionType type)
Select compression for sparse data types.
void SetPrintFactorStatistics(bool print_stat)
Set up verbose printing during the factor step.
void SetKrylovSolver(strumpack::KrylovSolver method)
Set the Krylov solver method to use.
void Mult(const Vector &x, Vector &y) const
Factor and solve the linear system .
void SetFromCommandLine()
Set options that were captured from the command line.
void SetPrintSolveStatistics(bool print_stat)
Set up verbose printing during the solve step.
void SetReorderingStrategy(strumpack::ReorderingStrategy method)
Set matrix reordering strategy.
Vector coefficient defined as a product of scalar and vector coefficients.
Scaled Operator B: x -> a A(x).
Definition operator.hpp:727
void SetColumnPermutation(superlu::ColPerm col_perm)
Specify how to permute the columns of the matrix.
Definition superlu.cpp:399
void Mult(const Vector &x, Vector &y) const
Factor and solve the linear system .
Definition superlu.cpp:587
void SetSymmetricPattern(bool sym)
Specify whether the matrix has a symmetric pattern to avoid extra work (default false)
Definition superlu.cpp:454
void SetOperator(const Operator &op)
Set the operator/matrix.
Definition superlu.cpp:475
void SetPrintStatistics(bool print_stat)
Specify whether to print the solver statistics (default true)
Definition superlu.cpp:385
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
prob_type
Definition ex25.cpp:149
real_t omega
Definition ex25p.cpp:141
void detJ_JT_J_inv_abs(const Vector &x, PML *pml, Vector &D)
Definition ex25p.cpp:947
void detJ_inv_JT_J_Im(const Vector &x, PML *pml, Vector &D)
Definition ex25p.cpp:989
T pow2(const T &x)
Definition ex25p.cpp:145
void E_bdr_data_Re(const Vector &x, Vector &E)
Definition ex25p.cpp:860
bool exact_known
Definition ex25p.cpp:143
void detJ_JT_J_inv_Re(const Vector &x, PML *pml, Vector &D)
Definition ex25p.cpp:913
int dim
Definition ex25p.cpp:142
real_t mu
Definition ex25p.cpp:139
real_t epsilon
Definition ex25p.cpp:140
void detJ_inv_JT_J_abs(const Vector &x, PML *pml, Vector &D)
Definition ex25p.cpp:1013
void maxwell_solution(const Vector &x, vector< complex< real_t > > &Eval)
Definition ex25p.cpp:738
void source(const Vector &x, Vector &f)
Definition ex25p.cpp:722
Array2D< real_t > comp_domain_bdr
Definition ex25p.cpp:136
void E_exact_Re(const Vector &x, Vector &E)
Definition ex25p.cpp:840
void detJ_JT_J_inv_Im(const Vector &x, PML *pml, Vector &D)
Definition ex25p.cpp:930
void detJ_inv_JT_J_Re(const Vector &x, PML *pml, Vector &D)
Definition ex25p.cpp:964
prob_type prob
Definition ex25p.cpp:155
prob_type
Definition ex25p.cpp:148
Definition ex25p.cpp:153
@ disc
Definition ex25p.cpp:150
@ lshape
Definition ex25p.cpp:151
@ beam
Definition ex25p.cpp:149
@ fichera
Definition ex25p.cpp:152
void E_exact_Im(const Vector &x, Vector &E)
Definition ex25p.cpp:850
Array2D< real_t > domain_bdr
Definition ex25p.cpp:137
void E_bdr_data_Im(const Vector &x, Vector &E)
Definition ex25p.cpp:887
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
@ PARMETIS
Sequential ordering on structure of using the PARMETIS package.
Definition superlu.hpp:71
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