MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pmaxwell.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11//
12// MFEM Ultraweak DPG Maxwell parallel example
13//
14// Compile with: make pmaxwell
15//
16// sample run
17// mpirun -np 4 pmaxwell -m ../../data/star.mesh -o 2 -sref 0 -pref 3 -rnum 0.5 -prob 0
18// mpirun -np 4 pmaxwell -m ../../data/inline-quad.mesh -o 3 -sref 0 -pref 3 -rnum 4.8 -sc -prob 0
19// mpirun -np 4 pmaxwell -m ../../data/inline-hex.mesh -o 2 -sref 0 -pref 1 -rnum 0.8 -sc -prob 0
20// mpirun -np 4 pmaxwell -m ../../data/inline-quad.mesh -o 3 -sref 1 -pref 3 -rnum 4.8 -sc -prob 2
21// mpirun -np 4 pmaxwell -o 3 -sref 1 -pref 2 -rnum 11.8 -sc -prob 3
22// mpirun -np 4 pmaxwell -o 3 -sref 1 -pref 2 -rnum 9.8 -sc -prob 4
23
24// AMR run. Note that this is a computationally intensive sample run.
25// We recommend trying it on a large machine with more mpi ranks
26// mpirun -np 4 pmaxwell -o 3 -sref 0 -pref 15 -prob 1 -theta 0.7 -sc
27
28// Description:
29// This example code demonstrates the use of MFEM to define and solve
30// the "ultraweak" (UW) DPG formulation for the Maxwell problem
31
32// ∇×(1/μ ∇×E) - ω² ϵ E = Ĵ , in Ω
33// E×n = E₀ , on ∂Ω
34
35// It solves the following kinds of problems
36// 1) Known exact solutions with error convergence rates
37// a) A manufactured solution problem where E is a plane beam
38// 2) Fichera "microwave" problem
39// 3) PML problems
40// a) Generic PML problem with point source given by the load
41// b) Plane wave scattering from a square
42// c) PML problem with a point source prescribed on the boundary
43
44// The DPG UW deals with the First Order System
45// i ω μ H + ∇ × E = 0, in Ω
46// -i ω ϵ E + ∇ × H = J, in Ω
47// E × n = E_0, on ∂Ω
48// Note: Ĵ = -iωJ
49
50// The ultraweak-DPG formulation is obtained by integration by parts of both
51// equations and the introduction of trace unknowns on the mesh skeleton
52
53// in 2D
54// E is vector valued and H is scalar.
55// (∇ × E, F) = (E, ∇ × F) + < n × E , F>
56// or (∇ ⋅ AE , F) = (AE, ∇ F) + < AE ⋅ n, F>
57// where A = [0 1; -1 0];
58
59// E ∈ (L²(Ω))² , H ∈ L²(Ω)
60// Ê ∈ H^-1/2(Ω)(Γₕ), Ĥ ∈ H^1/2(Γₕ)
61// i ω μ (H,F) + (E, ∇ × F) + < AÊ, F > = 0, ∀ F ∈ H¹
62// -i ω ϵ (E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω)
63// Ê = E₀ on ∂Ω
64// -------------------------------------------------------------------------
65// | | E | H | Ê | Ĥ | RHS |
66// -------------------------------------------------------------------------
67// | F | (E,∇ × F) | i ω μ (H,F) | < Ê, F > | | |
68// | | | | | | |
69// | G | -i ω ϵ (E,G) | (H,∇ × G) | | < Ĥ, G × n > | (J,G) |
70// where (F,G) ∈ H¹ × H(curl,Ω)
71
72// in 3D
73// E,H ∈ (L^2(Ω))³
74// Ê ∈ H_0^1/2(Ω)(curl, Γₕ), Ĥ ∈ H^-1/2(curl, Γₕ)
75// i ω μ (H,F) + (E,∇ × F) + < Ê, F × n > = 0, ∀ F ∈ H(curl,Ω)
76// -i ω ϵ (E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω)
77// Ê × n = E₀ on ∂Ω
78// -------------------------------------------------------------------------
79// | | E | H | Ê | Ĥ | RHS |
80// -------------------------------------------------------------------------
81// | F | (E,∇ × F) | i ω μ (H,F) | < n × Ê, F > | | |
82// | | | | | | |
83// | G | -i ω ϵ (E,G) | (H,∇ × G) | | < n × Ĥ, G > | (J,G) |
84// where (F,G) ∈ H(curl,Ω) × H(curl,Ω)
85
86// Here we use the "Adjoint Graph" norm on the test space i.e.,
87// ||(F,G)||²ᵥ = ||A^*(F,G)||² + ||(F,G)||² where A is the
88// maxwell operator defined by (1)
89
90// The PML formulation is
91
92// ∇×(1/μ α ∇×E) - ω² ϵ β E = Ĵ , in Ω
93// E×n = E₀ , on ∂Ω
94
95// where α = |J|⁻¹ Jᵀ J (in 2D it's the scalar |J|⁻¹),
96// β = |J| J⁻¹ J⁻ᵀ, J is the Jacobian of the stretching map
97// and |J| its determinant.
98
99// The first order system reads
100// i ω μ α⁻¹ H + ∇ × E = 0, in Ω
101// -i ω ϵ β E + ∇ × H = J, in Ω
102// E × n = E₀, on ∂Ω
103
104// and the ultraweak formulation is
105
106// in 2D
107// E ∈ (L²(Ω))² , H ∈ L²(Ω)
108// Ê ∈ H^-1/2(Ω)(Γₕ), Ĥ ∈ H^1/2(Γₕ)
109// i ω μ (α⁻¹ H,F) + (E, ∇ × F) + < AÊ, F > = 0, ∀ F ∈ H¹
110// -i ω ϵ (β E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω)
111// Ê = E₀ on ∂Ω
112// ---------------------------------------------------------------------------------
113// | | E | H | Ê | Ĥ | RHS |
114// ---------------------------------------------------------------------------------
115// | F | (E,∇ × F) | i ω μ (α⁻¹ H,F) | < Ê, F > | | |
116// | | | | | | |
117// | G | -i ω ϵ (β E,G) | (H,∇ × G) | | < Ĥ, G × n > | (J,G) |
118
119// where (F,G) ∈ H¹ × H(curl,Ω)
120
121//
122// in 3D
123// E,H ∈ (L^2(Ω))³
124// Ê ∈ H_0^1/2(Ω)(curl, Γ_h), Ĥ ∈ H^-1/2(curl, Γₕ)
125// i ω μ (α⁻¹ H,F) + (E,∇ × F) + < Ê, F × n > = 0, ∀ F ∈ H(curl,Ω)
126// -i ω ϵ (β E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω)
127// Ê × n = E_0 on ∂Ω
128// -------------------------------------------------------------------------------
129// | | E | H | Ê | Ĥ | RHS |
130// -------------------------------------------------------------------------------
131// | F | ( E,∇ × F) | i ω μ (α⁻¹ H,F) | < n × Ê, F > | | |
132// | | | | | | |
133// | G | -iωϵ (β E,G) | (H,∇ × G) | | < n × Ĥ, G > | (J,G) |
134// where (F,G) ∈ H(curl,Ω) × H(curl,Ω)
135
136// For more information see https://doi.org/10.1016/j.camwa.2021.01.017
137
138#include "mfem.hpp"
140#include "util/pml.hpp"
142#include <fstream>
143#include <iostream>
144
145using namespace std;
146using namespace mfem;
147using namespace mfem::common;
148
149void E_exact_r(const Vector &x, Vector & E_r);
150void E_exact_i(const Vector &x, Vector & E_i);
151
152void H_exact_r(const Vector &x, Vector & H_r);
153void H_exact_i(const Vector &x, Vector & H_i);
154
155
156void rhs_func_r(const Vector &x, Vector & J_r);
157void rhs_func_i(const Vector &x, Vector & J_i);
158
159void curlE_exact_r(const Vector &x, Vector &curlE_r);
160void curlE_exact_i(const Vector &x, Vector &curlE_i);
161void curlH_exact_r(const Vector &x,Vector &curlH_r);
162void curlH_exact_i(const Vector &x,Vector &curlH_i);
163
164void curlcurlE_exact_r(const Vector &x, Vector & curlcurlE_r);
165void curlcurlE_exact_i(const Vector &x, Vector & curlcurlE_i);
166
167void hatE_exact_r(const Vector & X, Vector & hatE_r);
168void hatE_exact_i(const Vector & X, Vector & hatE_i);
169
170void hatH_exact_r(const Vector & X, Vector & hatH_r);
171void hatH_exact_i(const Vector & X, Vector & hatH_i);
172
175
176void maxwell_solution(const Vector & X,
177 std::vector<complex<real_t>> &E);
178
179void maxwell_solution_curl(const Vector & X,
180 std::vector<complex<real_t>> &curlE);
181
182void maxwell_solution_curlcurl(const Vector & X,
183 std::vector<complex<real_t>> &curlcurlE);
184
185void source_function(const Vector &x, Vector & f);
186
187int dim;
190real_t mu = 1.0;
192
201
202static const char *enum_str[] =
203{
204 "plane_wave",
205 "fichera_oven",
206 "pml_general",
207 "pml_plane_wave_scatter",
208 "pml_pointsource"
209};
210
212
213int main(int argc, char *argv[])
214{
215 Mpi::Init();
216 int myid = Mpi::WorldRank();
217 Hypre::Init();
218
219 const char *mesh_file = "../../data/inline-quad.mesh";
220 int order = 1;
221 int delta_order = 1;
222 real_t rnum=1.0;
223 real_t theta = 0.0;
224 bool static_cond = false;
225 int iprob = 0;
226 int sr = 0;
227 int pr = 1;
228 bool exact_known = false;
229 bool with_pml = false;
230 bool visualization = true;
231 bool paraview = false;
232
233 OptionsParser args(argc, argv);
234 args.AddOption(&mesh_file, "-m", "--mesh",
235 "Mesh file to use.");
236 args.AddOption(&order, "-o", "--order",
237 "Finite element order (polynomial degree)");
238 args.AddOption(&rnum, "-rnum", "--number-of-wavelengths",
239 "Number of wavelengths");
240 args.AddOption(&mu, "-mu", "--permeability",
241 "Permeability of free space (or 1/(spring constant)).");
242 args.AddOption(&epsilon, "-eps", "--permittivity",
243 "Permittivity of free space (or mass constant).");
244 args.AddOption(&iprob, "-prob", "--problem", "Problem case"
245 " 0: plane wave, 1: Fichera 'oven', "
246 " 2: Generic PML problem with point source given as a load "
247 " 3: Scattering of a plane wave, "
248 " 4: Point source given on the boundary");
249 args.AddOption(&delta_order, "-do", "--delta-order",
250 "Order enrichment for DPG test space.");
251 args.AddOption(&theta, "-theta", "--theta",
252 "Theta parameter for AMR");
253 args.AddOption(&sr, "-sref", "--serial-ref",
254 "Number of parallel refinements.");
255 args.AddOption(&pr, "-pref", "--parallel-ref",
256 "Number of parallel refinements.");
257 args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
258 "--no-static-condensation", "Enable static condensation.");
259 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
260 "--no-visualization",
261 "Enable or disable GLVis visualization.");
262 args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
263 "--no-paraview",
264 "Enable or disable ParaView visualization.");
265 args.Parse();
266 if (!args.Good())
267 {
268 if (myid == 0)
269 {
270 args.PrintUsage(cout);
271 }
272 return 1;
273 }
274
275 if (iprob > 4) { iprob = 0; }
276 prob = (prob_type)iprob;
277 omega = 2.*M_PI*rnum;
278
279 if (prob == 0)
280 {
281 exact_known = true;
282 }
283 else if (prob == 1)
284 {
285 mesh_file = "meshes/fichera-waveguide.mesh";
286 omega = 5.0;
287 rnum = omega/(2.*M_PI);
288 }
289 else if (prob == 2)
290 {
291 with_pml = true;
292 }
293 else
294 {
295 with_pml = true;
296 mesh_file = "meshes/scatter.mesh";
297 }
298
299 if (myid == 0)
300 {
301 args.PrintOptions(cout);
302 }
303
304 Mesh mesh(mesh_file, 1, 1);
305 dim = mesh.Dimension();
306 MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
307
308 dimc = (dim == 3) ? 3 : 1;
309
310 for (int i = 0; i<sr; i++)
311 {
312 mesh.UniformRefinement();
313 }
314 mesh.EnsureNCMesh(false);
315
316 CartesianPML * pml = nullptr;
317 if (with_pml)
318 {
319 Array2D<real_t> length(dim, 2); length = 0.25;
320 pml = new CartesianPML(&mesh,length);
321 pml->SetOmega(omega);
323 }
324
325 ParMesh pmesh(MPI_COMM_WORLD, mesh);
326 mesh.Clear();
327
328 // PML element attribute marker
329 Array<int> attr;
330 Array<int> attrPML;
331 if (pml) { pml->SetAttributes(&pmesh, &attr, &attrPML); }
332
333 // Define spaces
334 enum TrialSpace
335 {
336 E_space = 0,
337 H_space = 1,
338 hatE_space = 2,
339 hatH_space = 3
340 };
341 enum TestSpace
342 {
343 F_space = 0,
344 G_space = 1
345 };
346 // L2 space for E
347 FiniteElementCollection *E_fec = new L2_FECollection(order-1,dim);
348 ParFiniteElementSpace *E_fes = new ParFiniteElementSpace(&pmesh,E_fec,dim);
349
350 // Vector L2 space for H
351 FiniteElementCollection *H_fec = new L2_FECollection(order-1,dim);
352 ParFiniteElementSpace *H_fes = new ParFiniteElementSpace(&pmesh,H_fec, dimc);
353
354 // H^-1/2 (curl) space for Ê
355 FiniteElementCollection * hatE_fec = nullptr;
356 FiniteElementCollection * hatH_fec = nullptr;
357 FiniteElementCollection * F_fec = nullptr;
358 int test_order = order+delta_order;
359 if (dim == 3)
360 {
361 hatE_fec = new ND_Trace_FECollection(order,dim);
362 hatH_fec = new ND_Trace_FECollection(order,dim);
363 F_fec = new ND_FECollection(test_order, dim);
364 }
365 else
366 {
367 hatE_fec = new RT_Trace_FECollection(order-1,dim);
368 hatH_fec = new H1_Trace_FECollection(order,dim);
369 F_fec = new H1_FECollection(test_order, dim);
370 }
371 ParFiniteElementSpace *hatE_fes = new ParFiniteElementSpace(&pmesh,hatE_fec);
372 ParFiniteElementSpace *hatH_fes = new ParFiniteElementSpace(&pmesh,hatH_fec);
373 FiniteElementCollection * G_fec = new ND_FECollection(test_order, dim);
374
377 trial_fes.Append(E_fes);
378 trial_fes.Append(H_fes);
379 trial_fes.Append(hatE_fes);
380 trial_fes.Append(hatH_fes);
381 test_fec.Append(F_fec);
382 test_fec.Append(G_fec);
383
384 // Bilinear form coefficients
385 ConstantCoefficient one(1.0);
389 ConstantCoefficient negepsomeg(-epsilon*omega);
391 ConstantCoefficient negmuomeg(-mu*omega);
392 // for the 2D case
393 DenseMatrix rot_mat(2);
394 rot_mat(0,0) = 0.; rot_mat(0,1) = 1.;
395 rot_mat(1,0) = -1.; rot_mat(1,1) = 0.;
396 MatrixConstantCoefficient rot(rot_mat);
397 ScalarMatrixProductCoefficient epsrot(epsomeg,rot);
398 ScalarMatrixProductCoefficient negepsrot(negepsomeg,rot);
399
400 Coefficient * epsomeg_cf = nullptr;
401 Coefficient * negepsomeg_cf = nullptr;
402 Coefficient * eps2omeg2_cf = nullptr;
403 Coefficient * muomeg_cf = nullptr;
404 Coefficient * negmuomeg_cf = nullptr;
405 Coefficient * mu2omeg2_cf = nullptr;
406 MatrixCoefficient *epsrot_cf = nullptr;
407 MatrixCoefficient *negepsrot_cf = nullptr;
408
409 if (pml)
410 {
411 epsomeg_cf = new RestrictedCoefficient(epsomeg,attr);
412 negepsomeg_cf = new RestrictedCoefficient(negepsomeg,attr);
413 eps2omeg2_cf = new RestrictedCoefficient(eps2omeg2,attr);
414 muomeg_cf = new RestrictedCoefficient(muomeg,attr);
415 negmuomeg_cf = new RestrictedCoefficient(negmuomeg,attr);
416 mu2omeg2_cf = new RestrictedCoefficient(mu2omeg2,attr);
417 epsrot_cf = new MatrixRestrictedCoefficient(epsrot,attr);
418 negepsrot_cf = new MatrixRestrictedCoefficient(negepsrot,attr);
419 }
420 else
421 {
422 epsomeg_cf = &epsomeg;
423 negepsomeg_cf = &negepsomeg;
424 eps2omeg2_cf = &eps2omeg2;
425 muomeg_cf = &muomeg;
426 negmuomeg_cf = &negmuomeg;
427 mu2omeg2_cf = &mu2omeg2;
428 epsrot_cf = &epsrot;
429 negepsrot_cf = &negepsrot;
430 }
431
432 // PML coefficients;
435 PmlCoefficient abs_detJ_2(abs_detJ_2_function,pml);
439
440 ProductCoefficient negmuomeg_detJ_r(negmuomeg,detJ_r);
441 ProductCoefficient negmuomeg_detJ_i(negmuomeg,detJ_i);
442 ProductCoefficient muomeg_detJ_r(muomeg,detJ_r);
443 ProductCoefficient mu2omeg2_detJ_2(mu2omeg2,abs_detJ_2);
444 ScalarMatrixProductCoefficient epsomeg_detJ_Jt_J_inv_i(epsomeg,
445 detJ_Jt_J_inv_i);
446 ScalarMatrixProductCoefficient epsomeg_detJ_Jt_J_inv_r(epsomeg,
447 detJ_Jt_J_inv_r);
448 ScalarMatrixProductCoefficient negepsomeg_detJ_Jt_J_inv_r(negepsomeg,
449 detJ_Jt_J_inv_r);
450 ScalarMatrixProductCoefficient muomeg_detJ_Jt_J_inv_r(muomeg,detJ_Jt_J_inv_r);
451 ScalarMatrixProductCoefficient negmuomeg_detJ_Jt_J_inv_i(negmuomeg,
452 detJ_Jt_J_inv_i);
453 ScalarMatrixProductCoefficient negmuomeg_detJ_Jt_J_inv_r(negmuomeg,
454 detJ_Jt_J_inv_r);
455 ScalarMatrixProductCoefficient mu2omeg2_detJ_Jt_J_inv_2(mu2omeg2,
456 abs_detJ_Jt_J_inv_2);
457 ScalarMatrixProductCoefficient eps2omeg2_detJ_Jt_J_inv_2(eps2omeg2,
458 abs_detJ_Jt_J_inv_2);
459
460 RestrictedCoefficient negmuomeg_detJ_r_restr(negmuomeg_detJ_r,attrPML);
461 RestrictedCoefficient negmuomeg_detJ_i_restr(negmuomeg_detJ_i,attrPML);
462 RestrictedCoefficient muomeg_detJ_r_restr(muomeg_detJ_r,attrPML);
463 RestrictedCoefficient mu2omeg2_detJ_2_restr(mu2omeg2_detJ_2,attrPML);
464 MatrixRestrictedCoefficient epsomeg_detJ_Jt_J_inv_i_restr(
465 epsomeg_detJ_Jt_J_inv_i,attrPML);
466 MatrixRestrictedCoefficient epsomeg_detJ_Jt_J_inv_r_restr(
467 epsomeg_detJ_Jt_J_inv_r,attrPML);
468 MatrixRestrictedCoefficient negepsomeg_detJ_Jt_J_inv_r_restr(
469 negepsomeg_detJ_Jt_J_inv_r,attrPML);
470 MatrixRestrictedCoefficient muomeg_detJ_Jt_J_inv_r_restr(muomeg_detJ_Jt_J_inv_r,
471 attrPML);
472 MatrixRestrictedCoefficient negmuomeg_detJ_Jt_J_inv_i_restr(
473 negmuomeg_detJ_Jt_J_inv_i,attrPML);
474 MatrixRestrictedCoefficient negmuomeg_detJ_Jt_J_inv_r_restr(
475 negmuomeg_detJ_Jt_J_inv_r,attrPML);
476 MatrixRestrictedCoefficient mu2omeg2_detJ_Jt_J_inv_2_restr(
477 mu2omeg2_detJ_Jt_J_inv_2,attrPML);
478 MatrixRestrictedCoefficient eps2omeg2_detJ_Jt_J_inv_2_restr(
479 eps2omeg2_detJ_Jt_J_inv_2,attrPML);
480
481 MatrixProductCoefficient * epsomeg_detJ_Jt_J_inv_i_rot = nullptr;
482 MatrixProductCoefficient * epsomeg_detJ_Jt_J_inv_r_rot = nullptr;
483 MatrixProductCoefficient * negepsomeg_detJ_Jt_J_inv_r_rot = nullptr;
484 MatrixRestrictedCoefficient * epsomeg_detJ_Jt_J_inv_i_rot_restr = nullptr;
485 MatrixRestrictedCoefficient * epsomeg_detJ_Jt_J_inv_r_rot_restr = nullptr;
486 MatrixRestrictedCoefficient * negepsomeg_detJ_Jt_J_inv_r_rot_restr = nullptr;
487
488 if (pml && dim == 2)
489 {
490 epsomeg_detJ_Jt_J_inv_i_rot = new MatrixProductCoefficient(
491 epsomeg_detJ_Jt_J_inv_i, rot);
492 epsomeg_detJ_Jt_J_inv_r_rot = new MatrixProductCoefficient(
493 epsomeg_detJ_Jt_J_inv_r, rot);
494 negepsomeg_detJ_Jt_J_inv_r_rot = new MatrixProductCoefficient(
495 negepsomeg_detJ_Jt_J_inv_r, rot);
496 epsomeg_detJ_Jt_J_inv_i_rot_restr = new MatrixRestrictedCoefficient(
497 *epsomeg_detJ_Jt_J_inv_i_rot, attrPML);
498 epsomeg_detJ_Jt_J_inv_r_rot_restr = new MatrixRestrictedCoefficient(
499 *epsomeg_detJ_Jt_J_inv_r_rot, attrPML);
500 negepsomeg_detJ_Jt_J_inv_r_rot_restr = new MatrixRestrictedCoefficient(
501 *negepsomeg_detJ_Jt_J_inv_r_rot, attrPML);
502 }
503
504 ParComplexDPGWeakForm * a = new ParComplexDPGWeakForm(trial_fes,test_fec);
505 a->StoreMatrices(); // needed for AMR
506
507 // (E,∇ × F)
508 a->AddTrialIntegrator(new TransposeIntegrator(new MixedCurlIntegrator(one)),
509 nullptr,TrialSpace::E_space, TestSpace::F_space);
510 // -i ω ϵ (E , G) = i (- ω ϵ E, G)
511 a->AddTrialIntegrator(nullptr,
512 new TransposeIntegrator(new VectorFEMassIntegrator(*negepsomeg_cf)),
513 TrialSpace::E_space,TestSpace::G_space);
514 // (H,∇ × G)
515 a->AddTrialIntegrator(new TransposeIntegrator(new MixedCurlIntegrator(one)),
516 nullptr,TrialSpace::H_space, TestSpace::G_space);
517 // < n×Ĥ ,G>
518 a->AddTrialIntegrator(new TangentTraceIntegrator,nullptr,
519 TrialSpace::hatH_space, TestSpace::G_space);
520 // test integrators
521 // (∇×G ,∇× δG)
522 a->AddTestIntegrator(new CurlCurlIntegrator(one),nullptr,
523 TestSpace::G_space,TestSpace::G_space);
524 // (G,δG)
525 a->AddTestIntegrator(new VectorFEMassIntegrator(one),nullptr,
526 TestSpace::G_space,TestSpace::G_space);
527
528 if (dim == 3)
529 {
530 // i ω μ (H, F)
531 a->AddTrialIntegrator(nullptr, new TransposeIntegrator(
532 new VectorFEMassIntegrator(*muomeg_cf)),
533 TrialSpace::H_space,TestSpace::F_space);
534 // < n×Ê,F>
535 a->AddTrialIntegrator(new TangentTraceIntegrator,nullptr,
536 TrialSpace::hatE_space, TestSpace::F_space);
537
538 // test integrators
539 // (∇×F,∇×δF)
540 a->AddTestIntegrator(new CurlCurlIntegrator(one),nullptr,
541 TestSpace::F_space, TestSpace::F_space);
542 // (F,δF)
543 a->AddTestIntegrator(new VectorFEMassIntegrator(one),nullptr,
544 TestSpace::F_space,TestSpace::F_space);
545 // μ^2 ω^2 (F,δF)
546 a->AddTestIntegrator(new VectorFEMassIntegrator(*mu2omeg2_cf),nullptr,
547 TestSpace::F_space, TestSpace::F_space);
548 // -i ω μ (F,∇ × δG) = i (F, -ω μ ∇ × δ G)
549 a->AddTestIntegrator(nullptr,new MixedVectorWeakCurlIntegrator(*negmuomeg_cf),
550 TestSpace::F_space, TestSpace::G_space);
551 // -i ω ϵ (∇ × F, δG)
552 a->AddTestIntegrator(nullptr,new MixedVectorCurlIntegrator(*negepsomeg_cf),
553 TestSpace::F_space, TestSpace::G_space);
554 // i ω μ (∇ × G,δF)
555 a->AddTestIntegrator(nullptr,new MixedVectorCurlIntegrator(*muomeg_cf),
556 TestSpace::G_space, TestSpace::F_space);
557 // i ω ϵ (G, ∇ × δF )
558 a->AddTestIntegrator(nullptr,new MixedVectorWeakCurlIntegrator(*epsomeg_cf),
559 TestSpace::G_space, TestSpace::F_space);
560 // ϵ^2 ω^2 (G,δG)
561 a->AddTestIntegrator(new VectorFEMassIntegrator(*eps2omeg2_cf),nullptr,
562 TestSpace::G_space, TestSpace::G_space);
563 }
564 else
565 {
566 // i ω μ (H, F)
567 a->AddTrialIntegrator(nullptr,new MixedScalarMassIntegrator(*muomeg_cf),
568 TrialSpace::H_space, TestSpace::F_space);
569 // < n×Ê,F>
570 a->AddTrialIntegrator(new TraceIntegrator,nullptr,
571 TrialSpace::hatE_space, TestSpace::F_space);
572 // test integrators
573 // (∇F,∇δF)
574 a->AddTestIntegrator(new DiffusionIntegrator(one),nullptr,
575 TestSpace::F_space, TestSpace::F_space);
576 // (F,δF)
577 a->AddTestIntegrator(new MassIntegrator(one),nullptr,
578 TestSpace::F_space, TestSpace::F_space);
579 // μ^2 ω^2 (F,δF)
580 a->AddTestIntegrator(new MassIntegrator(*mu2omeg2_cf),nullptr,
581 TestSpace::F_space, TestSpace::F_space);
582 // -i ω μ (F,∇ × δG) = i (F, -ω μ ∇ × δ G)
583 a->AddTestIntegrator(nullptr,
584 new TransposeIntegrator(new MixedCurlIntegrator(*negmuomeg_cf)),
585 TestSpace::F_space, TestSpace::G_space);
586 // -i ω ϵ (∇ × F, δG) = i (- ω ϵ A ∇ F,δG), A = [0 1; -1; 0]
587 a->AddTestIntegrator(nullptr,new MixedVectorGradientIntegrator(*negepsrot_cf),
588 TestSpace::F_space, TestSpace::G_space);
589 // i ω μ (∇ × G,δF) = i (ω μ ∇ × G, δF )
590 a->AddTestIntegrator(nullptr,new MixedCurlIntegrator(*muomeg_cf),
591 TestSpace::G_space, TestSpace::F_space);
592 // i ω ϵ (G, ∇ × δF ) = i (ω ϵ G, A ∇ δF) = i ( G , ω ϵ A ∇ δF)
593 a->AddTestIntegrator(nullptr,
595 new MixedVectorGradientIntegrator(*epsrot_cf)),
596 TestSpace::G_space, TestSpace::F_space);
597 // ϵ^2 ω^2 (G, δG)
598 a->AddTestIntegrator(new VectorFEMassIntegrator(*eps2omeg2_cf),nullptr,
599 TestSpace::G_space, TestSpace::G_space);
600 }
601 if (pml)
602 {
603 //trial integrators
604 // -i ω ϵ (β E , G) = -i ω ϵ ((β_re + i β_im) E, G)
605 // = (ω ϵ β_im E, G) + i (- ω ϵ β_re E, G)
606 a->AddTrialIntegrator(
608 epsomeg_detJ_Jt_J_inv_i_restr)),
610 negepsomeg_detJ_Jt_J_inv_r_restr)),
611 TrialSpace::E_space,TestSpace::G_space);
612 if (dim == 3)
613 {
614 //trial integrators
615 // i ω μ (α^-1 H, F) = i ω μ ( (α^-1_re + i α^-1_im) H, F)
616 // = (- ω μ α^-1_im, H,F) + i *(ω μ α^-1_re H, F)
617 a->AddTrialIntegrator(
619 negmuomeg_detJ_Jt_J_inv_i_restr)),
621 muomeg_detJ_Jt_J_inv_r_restr)),
622 TrialSpace::H_space, TestSpace::F_space);
623 // test integrators
624 // μ^2 ω^2 (|α|^-2 F,δF)
625 a->AddTestIntegrator(
626 new VectorFEMassIntegrator(mu2omeg2_detJ_Jt_J_inv_2_restr),nullptr,
627 TestSpace::F_space, TestSpace::F_space);
628 // -i ω μ (α^-* F,∇ × δG) = i (F, - ω μ α^-1 ∇ × δ G)
629 // = i (F, - ω μ (α^-1_re + i α^-1_im) ∇ × δ G)
630 // = (F, - ω μ α^-1_im ∇ × δ G) + i (F, - ω μ α^-1_re ∇×δG)
631 a->AddTestIntegrator(new MixedVectorWeakCurlIntegrator(
632 negmuomeg_detJ_Jt_J_inv_i_restr),
633 new MixedVectorWeakCurlIntegrator(negmuomeg_detJ_Jt_J_inv_r_restr),
634 TestSpace::F_space,TestSpace::G_space);
635 // -i ω ϵ (β ∇ × F, δG) = -i ω ϵ ((β_re + i β_im) ∇ × F, δG)
636 // = (ω ϵ β_im ∇ × F, δG) + i (- ω ϵ β_re ∇ × F, δG)
637 a->AddTestIntegrator(new MixedVectorCurlIntegrator(
638 epsomeg_detJ_Jt_J_inv_i_restr),
639 new MixedVectorCurlIntegrator(negepsomeg_detJ_Jt_J_inv_r_restr),
640 TestSpace::F_space,TestSpace::G_space);
641 // i ω μ (α^-1 ∇ × G,δF) = i ω μ ((α^-1_re + i α^-1_im) ∇ × G,δF)
642 // = (- ω μ α^-1_im ∇ × G,δF) + i (ω μ α^-1_re ∇ × G,δF)
643 a->AddTestIntegrator(new MixedVectorCurlIntegrator(
644 negmuomeg_detJ_Jt_J_inv_i_restr),
645 new MixedVectorCurlIntegrator(muomeg_detJ_Jt_J_inv_r_restr),
646 TestSpace::G_space, TestSpace::F_space);
647 // i ω ϵ (β^* G, ∇×δF) = i ω ϵ ( (β_re - i β_im) G, ∇×δF)
648 // = (ω ϵ β_im G, ∇×δF) + i ( ω ϵ β_re G, ∇×δF)
649 a->AddTestIntegrator(new MixedVectorWeakCurlIntegrator(
650 epsomeg_detJ_Jt_J_inv_i_restr),
651 new MixedVectorWeakCurlIntegrator(epsomeg_detJ_Jt_J_inv_r_restr),
652 TestSpace::G_space, TestSpace::F_space);
653 // ϵ^2 ω^2 (|β|^2 G,δG)
654 a->AddTestIntegrator(new VectorFEMassIntegrator(
655 eps2omeg2_detJ_Jt_J_inv_2_restr),nullptr,
656 TestSpace::G_space, TestSpace::G_space);
657 }
658 else
659 {
660 //trial integrators
661 // i ω μ (α^-1 H, F) = i ω μ ( (α^-1_re + i α^-1_im) H, F)
662 // = (- ω μ α^-1_im, H,F) + i *(ω μ α^-1_re H, F)
663 a->AddTrialIntegrator(
664 new MixedScalarMassIntegrator(negmuomeg_detJ_i_restr),
665 new MixedScalarMassIntegrator(muomeg_detJ_r_restr),
666 TrialSpace::H_space, TestSpace::F_space);
667 // test integrators
668 // μ^2 ω^2 (|α|^-2 F,δF)
669 a->AddTestIntegrator(new MassIntegrator(mu2omeg2_detJ_2_restr),nullptr,
670 TestSpace::F_space, TestSpace::F_space);
671 // -i ω μ (α^-* F,∇ × δG) = (F, ω μ α^-1 ∇ × δ G)
672 // =(F, - ω μ α^-1_im ∇ × δ G) + i (F, - ω μ α^-1_re ∇×δG)
673 a->AddTestIntegrator(
674 new TransposeIntegrator(new MixedCurlIntegrator(negmuomeg_detJ_i_restr)),
675 new TransposeIntegrator(new MixedCurlIntegrator(negmuomeg_detJ_r_restr)),
676 TestSpace::F_space, TestSpace::G_space);
677 // -i ω ϵ (β ∇ × F, δG) = i (- ω ϵ β A ∇ F,δG), A = [0 1; -1; 0]
678 // = (ω ϵ β_im A ∇ F, δG) + i (- ω ϵ β_re A ∇ F, δG)
679 a->AddTestIntegrator(new MixedVectorGradientIntegrator(
680 *epsomeg_detJ_Jt_J_inv_i_rot_restr),
681 new MixedVectorGradientIntegrator(*negepsomeg_detJ_Jt_J_inv_r_rot_restr),
682 TestSpace::F_space, TestSpace::G_space);
683 // i ω μ (α^-1 ∇ × G,δF) = i (ω μ α^-1 ∇ × G, δF )
684 // = (- ω μ α^-1_im ∇ × G,δF) + i (ω μ α^-1_re ∇ × G,δF)
685 a->AddTestIntegrator(new MixedCurlIntegrator(negmuomeg_detJ_i_restr),
686 new MixedCurlIntegrator(muomeg_detJ_r_restr),
687 TestSpace::G_space, TestSpace::F_space);
688 // i ω ϵ (β^* G, ∇ × δF ) = i ( G , ω ϵ β A ∇ δF)
689 // = ( G , ω ϵ β_im A ∇ δF) + i ( G , ω ϵ β_re A ∇ δF)
690 a->AddTestIntegrator(
692 *epsomeg_detJ_Jt_J_inv_i_rot_restr)),
694 *epsomeg_detJ_Jt_J_inv_r_rot_restr)),
695 TestSpace::G_space, TestSpace::F_space);
696 // ϵ^2 ω^2 (|β|^2 G,δG)
697 a->AddTestIntegrator(new VectorFEMassIntegrator(
698 eps2omeg2_detJ_Jt_J_inv_2_restr),nullptr,
699 TestSpace::G_space, TestSpace::G_space);
700 }
701 }
702 // RHS
706 if (prob == 0)
707 {
708 a->AddDomainLFIntegrator(new VectorFEDomainLFIntegrator(f_rhs_r),
709 new VectorFEDomainLFIntegrator(f_rhs_i),
710 TestSpace::G_space);
711 }
712 else if (prob == 2)
713 {
714 a->AddDomainLFIntegrator(new VectorFEDomainLFIntegrator(f_source),nullptr,
715 TestSpace::G_space);
716 }
717
720
721 socketstream E_out_r;
722 socketstream H_out_r;
723 if (myid == 0)
724 {
725 std::cout << "\n Ref |"
726 << " Dofs |"
727 << " ω |" ;
728 if (exact_known)
729 {
730 std::cout << " L2 Error |"
731 << " Rate |" ;
732 }
733 std::cout << " Residual |"
734 << " Rate |"
735 << " PCG it |" << endl;
736 std::cout << std::string((exact_known) ? 82 : 60,'-')
737 << endl;
738 }
739
740 real_t res0 = 0.;
741 real_t err0 = 0.;
742 int dof0 = 0; // init to suppress gcc warning
743
744 Array<int> elements_to_refine;
745
746 ParGridFunction E_r, E_i, H_r, H_i;
747
748 ParaViewDataCollection * paraview_dc = nullptr;
749
750 if (paraview)
751 {
752 paraview_dc = new ParaViewDataCollection(enum_str[prob], &pmesh);
753 paraview_dc->SetPrefixPath("ParaView/Maxwell");
754 paraview_dc->SetLevelsOfDetail(order);
755 paraview_dc->SetCycle(0);
756 paraview_dc->SetDataFormat(VTKFormat::BINARY);
757 paraview_dc->SetHighOrderOutput(true);
758 paraview_dc->SetTime(0.0); // set the time
759 paraview_dc->RegisterField("E_r",&E_r);
760 paraview_dc->RegisterField("E_i",&E_i);
761 paraview_dc->RegisterField("H_r",&H_r);
762 paraview_dc->RegisterField("H_i",&H_i);
763 }
764
765 if (static_cond) { a->EnableStaticCondensation(); }
766 for (int it = 0; it<=pr; it++)
767 {
768 a->Assemble();
769
770 Array<int> ess_tdof_list;
771 Array<int> ess_bdr;
772 if (pmesh.bdr_attributes.Size())
773 {
774 ess_bdr.SetSize(pmesh.bdr_attributes.Max());
775 ess_bdr = 1;
776 hatE_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
777 if (pml)
778 {
779 ess_bdr = 0;
780 ess_bdr[1] = 1;
781 }
782 }
783
784 // shift the ess_tdofs
785 for (int j = 0; j < ess_tdof_list.Size(); j++)
786 {
787 ess_tdof_list[j] += E_fes->GetTrueVSize() + H_fes->GetTrueVSize();
788 }
789
790 Array<int> offsets(5);
791 offsets[0] = 0;
792 offsets[1] = E_fes->GetVSize();
793 offsets[2] = H_fes->GetVSize();
794 offsets[3] = hatE_fes->GetVSize();
795 offsets[4] = hatH_fes->GetVSize();
796 offsets.PartialSum();
797
798 Vector x(2*offsets.Last());
799 x = 0.;
800
801 if (prob != 2)
802 {
803 ParGridFunction hatE_gf_r(hatE_fes, x, offsets[2]);
804 ParGridFunction hatE_gf_i(hatE_fes, x, offsets.Last() + offsets[2]);
805 if (dim == 3)
806 {
807 hatE_gf_r.ProjectBdrCoefficientTangent(hatEex_r, ess_bdr);
808 hatE_gf_i.ProjectBdrCoefficientTangent(hatEex_i, ess_bdr);
809 }
810 else
811 {
812 hatE_gf_r.ProjectBdrCoefficientNormal(hatEex_r, ess_bdr);
813 hatE_gf_i.ProjectBdrCoefficientNormal(hatEex_i, ess_bdr);
814 }
815 }
816
817 OperatorPtr Ah;
818 Vector X,B;
819 a->FormLinearSystem(ess_tdof_list,x,Ah, X,B);
820
821 ComplexOperator * Ahc = Ah.As<ComplexOperator>();
822
823 BlockOperator * BlockA_r = dynamic_cast<BlockOperator *>(&Ahc->real());
824 BlockOperator * BlockA_i = dynamic_cast<BlockOperator *>(&Ahc->imag());
825
826 int num_blocks = BlockA_r->NumRowBlocks();
827 Array<int> tdof_offsets(2*num_blocks+1);
828
829 tdof_offsets[0] = 0;
830 int skip = (static_cond) ? 0 : 2;
831 int k = (static_cond) ? 2 : 0;
832 for (int i=0; i<num_blocks; i++)
833 {
834 tdof_offsets[i+1] = trial_fes[i+k]->GetTrueVSize();
835 tdof_offsets[num_blocks+i+1] = trial_fes[i+k]->GetTrueVSize();
836 }
837 tdof_offsets.PartialSum();
838
839 BlockOperator blockA(tdof_offsets);
840 for (int i = 0; i<num_blocks; i++)
841 {
842 for (int j = 0; j<num_blocks; j++)
843 {
844 blockA.SetBlock(i,j,&BlockA_r->GetBlock(i,j));
845 blockA.SetBlock(i,j+num_blocks,&BlockA_i->GetBlock(i,j), -1.0);
846 blockA.SetBlock(i+num_blocks,j+num_blocks,&BlockA_r->GetBlock(i,j));
847 blockA.SetBlock(i+num_blocks,j,&BlockA_i->GetBlock(i,j));
848 }
849 }
850
851 X = 0.;
852 BlockDiagonalPreconditioner M(tdof_offsets);
853
854 if (!static_cond)
855 {
856 HypreBoomerAMG * solver_E = new HypreBoomerAMG((HypreParMatrix &)
857 BlockA_r->GetBlock(0,0));
858 solver_E->SetPrintLevel(0);
859 solver_E->SetSystemsOptions(dim);
860 HypreBoomerAMG * solver_H = new HypreBoomerAMG((HypreParMatrix &)
861 BlockA_r->GetBlock(1,1));
862 solver_H->SetPrintLevel(0);
863 solver_H->SetSystemsOptions(dim);
864 M.SetDiagonalBlock(0,solver_E);
865 M.SetDiagonalBlock(1,solver_H);
866 M.SetDiagonalBlock(num_blocks,solver_E);
867 M.SetDiagonalBlock(num_blocks+1,solver_H);
868 }
869
870 HypreSolver * solver_hatH = nullptr;
871 HypreAMS * solver_hatE = new HypreAMS((HypreParMatrix &)BlockA_r->GetBlock(skip,
872 skip),
873 hatE_fes);
874 solver_hatE->SetPrintLevel(0);
875 if (dim == 2)
876 {
877 solver_hatH = new HypreBoomerAMG((HypreParMatrix &)BlockA_r->GetBlock(skip+1,
878 skip+1));
879 dynamic_cast<HypreBoomerAMG*>(solver_hatH)->SetPrintLevel(0);
880 }
881 else
882 {
883 solver_hatH = new HypreAMS((HypreParMatrix &)BlockA_r->GetBlock(skip+1,skip+1),
884 hatH_fes);
885 dynamic_cast<HypreAMS*>(solver_hatH)->SetPrintLevel(0);
886 }
887
888 M.SetDiagonalBlock(skip,solver_hatE);
889 M.SetDiagonalBlock(skip+1,solver_hatH);
890 M.SetDiagonalBlock(skip+num_blocks,solver_hatE);
891 M.SetDiagonalBlock(skip+num_blocks+1,solver_hatH);
892
893 CGSolver cg(MPI_COMM_WORLD);
894 cg.SetRelTol(1e-6);
895 cg.SetMaxIter(10000);
896 cg.SetPrintLevel(0);
897 cg.SetPreconditioner(M);
898 cg.SetOperator(blockA);
899 cg.Mult(B, X);
900
901 for (int i = 0; i<num_blocks; i++)
902 {
903 delete &M.GetDiagonalBlock(i);
904 }
905
906 int num_iter = cg.GetNumIterations();
907
908 a->RecoverFEMSolution(X,x);
909
910 Vector & residuals = a->ComputeResidual(x);
911
912 real_t residual = residuals.Norml2();
913 real_t maxresidual = residuals.Max();
914 real_t globalresidual = residual * residual;
915 MPI_Allreduce(MPI_IN_PLACE, &maxresidual, 1, MPITypeMap<real_t>::mpi_type,
916 MPI_MAX, MPI_COMM_WORLD);
917 MPI_Allreduce(MPI_IN_PLACE, &globalresidual, 1,
918 MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
919
920 globalresidual = sqrt(globalresidual);
921
922 E_r.MakeRef(E_fes,x, 0);
923 E_i.MakeRef(E_fes,x, offsets.Last());
924
925 H_r.MakeRef(H_fes,x, offsets[1]);
926 H_i.MakeRef(H_fes,x, offsets.Last()+offsets[1]);
927
928 int dofs = 0;
929 for (int i = 0; i<trial_fes.Size(); i++)
930 {
931 dofs += trial_fes[i]->GlobalTrueVSize();
932 }
933
934 real_t L2Error = 0.0;
935 real_t rate_err = 0.0;
936
937 if (exact_known)
938 {
943 real_t E_err_r = E_r.ComputeL2Error(E_ex_r);
944 real_t E_err_i = E_i.ComputeL2Error(E_ex_i);
945 real_t H_err_r = H_r.ComputeL2Error(H_ex_r);
946 real_t H_err_i = H_i.ComputeL2Error(H_ex_i);
947 L2Error = sqrt( E_err_r*E_err_r + E_err_i*E_err_i
948 + H_err_r*H_err_r + H_err_i*H_err_i );
949 rate_err = (it) ? dim*log(err0/L2Error)/log((real_t)dof0/dofs) : 0.0;
950 err0 = L2Error;
951 }
952
953 real_t rate_res = (it) ? dim*log(res0/globalresidual)/log((
954 real_t)dof0/dofs) : 0.0;
955
956 res0 = globalresidual;
957 dof0 = dofs;
958
959 if (myid == 0)
960 {
961 std::ios oldState(nullptr);
962 oldState.copyfmt(std::cout);
963 std::cout << std::right << std::setw(5) << it << " | "
964 << std::setw(10) << dof0 << " | "
965 << std::setprecision(1) << std::fixed
966 << std::setw(4) << 2.0*rnum << " π | "
967 << std::setprecision(3);
968 if (exact_known)
969 {
970 std::cout << std::setw(10) << std::scientific << err0 << " | "
971 << std::setprecision(2)
972 << std::setw(6) << std::fixed << rate_err << " | " ;
973 }
974 std::cout << std::setprecision(3)
975 << std::setw(10) << std::scientific << res0 << " | "
976 << std::setprecision(2)
977 << std::setw(6) << std::fixed << rate_res << " | "
978 << std::setw(6) << std::fixed << num_iter << " | "
979 << std::endl;
980 std::cout.copyfmt(oldState);
981 }
982
983 if (visualization)
984 {
985 const char * keys = (it == 0 && dim == 2) ? "jRcml\n" : nullptr;
986 char vishost[] = "localhost";
987 int visport = 19916;
988 VisualizeField(E_out_r,vishost, visport, E_r,
989 "Numerical Electric field (real part)", 0, 0, 500, 500, keys);
990 VisualizeField(H_out_r,vishost, visport, H_r,
991 "Numerical Magnetic field (real part)", 501, 0, 500, 500, keys);
992 }
993
994 if (paraview)
995 {
996 paraview_dc->SetCycle(it);
997 paraview_dc->SetTime((real_t)it);
998 paraview_dc->Save();
999 }
1000
1001 if (it == pr)
1002 {
1003 break;
1004 }
1005
1006 if (theta > 0.0)
1007 {
1008 elements_to_refine.SetSize(0);
1009 for (int iel = 0; iel<pmesh.GetNE(); iel++)
1010 {
1011 if (residuals[iel] > theta * maxresidual)
1012 {
1013 elements_to_refine.Append(iel);
1014 }
1015 }
1016 pmesh.GeneralRefinement(elements_to_refine,1,1);
1017 }
1018 else
1019 {
1020 pmesh.UniformRefinement();
1021 }
1022 if (pml) { pml->SetAttributes(&pmesh); }
1023 for (int i =0; i<trial_fes.Size(); i++)
1024 {
1025 trial_fes[i]->Update(false);
1026 }
1027 a->Update();
1028 }
1029
1030 if (pml && dim == 2)
1031 {
1032 delete epsomeg_detJ_Jt_J_inv_i_rot;
1033 delete epsomeg_detJ_Jt_J_inv_r_rot;
1034 delete negepsomeg_detJ_Jt_J_inv_r_rot;
1035 delete epsomeg_detJ_Jt_J_inv_i_rot_restr;
1036 delete epsomeg_detJ_Jt_J_inv_r_rot_restr;
1037 delete negepsomeg_detJ_Jt_J_inv_r_rot_restr;
1038 }
1039
1040 if (paraview)
1041 {
1042 delete paraview_dc;
1043 }
1044
1045 delete a;
1046 delete F_fec;
1047 delete G_fec;
1048 delete hatH_fes;
1049 delete hatH_fec;
1050 delete hatE_fes;
1051 delete hatE_fec;
1052 delete H_fec;
1053 delete E_fec;
1054 delete H_fes;
1055 delete E_fes;
1056
1057 return 0;
1058}
1059
1060void E_exact_r(const Vector &x, Vector & E_r)
1061{
1062 std::vector<std::complex<real_t>> E;
1063 maxwell_solution(x,E);
1064 E_r.SetSize(E.size());
1065 for (unsigned i = 0; i < E.size(); i++)
1066 {
1067 E_r[i]= E[i].real();
1068 }
1069}
1070
1071void E_exact_i(const Vector &x, Vector & E_i)
1072{
1073 std::vector<std::complex<real_t>> E;
1074 maxwell_solution(x, E);
1075 E_i.SetSize(E.size());
1076 for (unsigned i = 0; i < E.size(); i++)
1077 {
1078 E_i[i]= E[i].imag();
1079 }
1080}
1081
1082void curlE_exact_r(const Vector &x, Vector &curlE_r)
1083{
1084 std::vector<std::complex<real_t>> curlE;
1085 maxwell_solution_curl(x, curlE);
1086 curlE_r.SetSize(curlE.size());
1087 for (unsigned i = 0; i < curlE.size(); i++)
1088 {
1089 curlE_r[i]= curlE[i].real();
1090 }
1091}
1092
1093void curlE_exact_i(const Vector &x, Vector &curlE_i)
1094{
1095 std::vector<std::complex<real_t>> curlE;
1096 maxwell_solution_curl(x, curlE);
1097 curlE_i.SetSize(curlE.size());
1098 for (unsigned i = 0; i < curlE.size(); i++)
1099 {
1100 curlE_i[i]= curlE[i].imag();
1101 }
1102}
1103
1104void curlcurlE_exact_r(const Vector &x, Vector & curlcurlE_r)
1105{
1106 std::vector<std::complex<real_t>> curlcurlE;
1107 maxwell_solution_curlcurl(x, curlcurlE);
1108 curlcurlE_r.SetSize(curlcurlE.size());
1109 for (unsigned i = 0; i < curlcurlE.size(); i++)
1110 {
1111 curlcurlE_r[i]= curlcurlE[i].real();
1112 }
1113}
1114
1115void curlcurlE_exact_i(const Vector &x, Vector & curlcurlE_i)
1116{
1117 std::vector<std::complex<real_t>> curlcurlE;
1118 maxwell_solution_curlcurl(x, curlcurlE);
1119 curlcurlE_i.SetSize(curlcurlE.size());
1120 for (unsigned i = 0; i < curlcurlE.size(); i++)
1121 {
1122 curlcurlE_i[i]= curlcurlE[i].imag();
1123 }
1124}
1125
1126
1127void H_exact_r(const Vector &x, Vector & H_r)
1128{
1129 // H = i ∇ × E / ω μ
1130 // H_r = - ∇ × E_i / ω μ
1131 Vector curlE_i;
1132 curlE_exact_i(x,curlE_i);
1133 H_r.SetSize(dimc);
1134 for (int i = 0; i<dimc; i++)
1135 {
1136 H_r(i) = - curlE_i(i) / (omega * mu);
1137 }
1138}
1139
1140void H_exact_i(const Vector &x, Vector & H_i)
1141{
1142 // H = i ∇ × E / ω μ
1143 // H_i = ∇ × E_r / ω μ
1144 Vector curlE_r;
1145 curlE_exact_r(x,curlE_r);
1146 H_i.SetSize(dimc);
1147 for (int i = 0; i<dimc; i++)
1148 {
1149 H_i(i) = curlE_r(i) / (omega * mu);
1150 }
1151}
1152
1153void curlH_exact_r(const Vector &x,Vector &curlH_r)
1154{
1155 // ∇ × H_r = - ∇ × ∇ × E_i / ω μ
1156 Vector curlcurlE_i;
1157 curlcurlE_exact_i(x,curlcurlE_i);
1158 curlH_r.SetSize(dim);
1159 for (int i = 0; i<dim; i++)
1160 {
1161 curlH_r(i) = -curlcurlE_i(i) / (omega * mu);
1162 }
1163}
1164
1165void curlH_exact_i(const Vector &x,Vector &curlH_i)
1166{
1167 // ∇ × H_i = ∇ × ∇ × E_r / ω μ
1168 Vector curlcurlE_r;
1169 curlcurlE_exact_r(x,curlcurlE_r);
1170 curlH_i.SetSize(dim);
1171 for (int i = 0; i<dim; i++)
1172 {
1173 curlH_i(i) = curlcurlE_r(i) / (omega * mu);
1174 }
1175}
1176
1177void hatE_exact_r(const Vector & x, Vector & hatE_r)
1178{
1179 if (dim == 3)
1180 {
1181 E_exact_r(x,hatE_r);
1182 }
1183 else
1184 {
1185 Vector E_r;
1186 E_exact_r(x,E_r);
1187 hatE_r.SetSize(hatE_r.Size());
1188 // rotate E_hat
1189 hatE_r[0] = E_r[1];
1190 hatE_r[1] = -E_r[0];
1191 }
1192}
1193
1194void hatE_exact_i(const Vector & x, Vector & hatE_i)
1195{
1196 if (dim == 3)
1197 {
1198 E_exact_i(x,hatE_i);
1199 }
1200 else
1201 {
1202 Vector E_i;
1203 E_exact_i(x,E_i);
1204 hatE_i.SetSize(hatE_i.Size());
1205 // rotate E_hat
1206 hatE_i[0] = E_i[1];
1207 hatE_i[1] = -E_i[0];
1208 }
1209}
1210
1211void hatH_exact_r(const Vector & x, Vector & hatH_r)
1212{
1213 H_exact_r(x,hatH_r);
1214}
1215
1216void hatH_exact_i(const Vector & x, Vector & hatH_i)
1217{
1218 H_exact_i(x,hatH_i);
1219}
1220
1222{
1223 Vector hatH_r;
1224 H_exact_r(x,hatH_r);
1225 return hatH_r[0];
1226}
1227
1229{
1230 Vector hatH_i;
1231 H_exact_i(x,hatH_i);
1232 return hatH_i[0];
1233}
1234
1235// J = -i ω ϵ E + ∇ × H
1236// J_r + iJ_i = -i ω ϵ (E_r + i E_i) + ∇ × (H_r + i H_i)
1237void rhs_func_r(const Vector &x, Vector & J_r)
1238{
1239 // J_r = ω ϵ E_i + ∇ × H_r
1240 Vector E_i, curlH_r;
1241 E_exact_i(x,E_i);
1242 curlH_exact_r(x,curlH_r);
1243 J_r.SetSize(dim);
1244 for (int i = 0; i<dim; i++)
1245 {
1246 J_r(i) = omega * epsilon * E_i(i) + curlH_r(i);
1247 }
1248}
1249
1250void rhs_func_i(const Vector &x, Vector & J_i)
1251{
1252 // J_i = - ω ϵ E_r + ∇ × H_i
1253 Vector E_r, curlH_i;
1254 E_exact_r(x,E_r);
1255 curlH_exact_i(x,curlH_i);
1256 J_i.SetSize(dim);
1257 for (int i = 0; i<dim; i++)
1258 {
1259 J_i(i) = -omega * epsilon * E_r(i) + curlH_i(i);
1260 }
1261}
1262
1263void maxwell_solution(const Vector & X, std::vector<complex<real_t>> &E)
1264{
1265 complex<real_t> zi = complex<real_t>(0., 1.);
1266 E.resize(dim);
1267 for (int i = 0; i < dim; ++i)
1268 {
1269 E[i] = 0.0;
1270 }
1271 switch (prob)
1272 {
1273 case plane_wave:
1274 {
1275 E[0] = exp(zi * omega * (X.Sum()));
1276 }
1277 break;
1279 {
1280 E[1] = exp(zi * omega * (X(0)));
1281 }
1282 break;
1283 case fichera_oven:
1284 {
1285 if (abs(X(2) - 3.0) < 1e-10)
1286 {
1287 E[0] = sin(M_PI*X(1));
1288 }
1289 }
1290 break;
1291 case pml_pointsource:
1292 {
1293 Vector shift(dim);
1294 real_t k = omega * sqrt(epsilon * mu);
1295 shift = -0.5;
1296
1297 if (dim == 2)
1298 {
1299 real_t x0 = X(0) + shift(0);
1300 real_t x1 = X(1) + shift(1);
1301 real_t r = sqrt(x0 * x0 + x1 * x1);
1302 real_t beta = k * r;
1303
1304 // Bessel functions
1305 complex<real_t> Ho, Ho_r, Ho_rr;
1306 Ho = real_t(jn(0, beta)) + zi * real_t(yn(0, beta));
1307 Ho_r = -k * real_t(jn(1, beta)) + zi * real_t(yn(1, beta));
1308 Ho_rr = -k * k * (1_r / beta *
1309 (real_t(jn(1, beta)) + zi * real_t(yn(1, beta))) -
1310 (real_t(jn(2, beta)) + zi * real_t(yn(2, beta))));
1311
1312 // First derivatives
1313 real_t r_x = x0 / r;
1314 real_t r_y = x1 / r;
1315 real_t r_xy = -(r_x / r) * r_y;
1316 real_t r_xx = (1.0 / r) * (1.0 - r_x * r_x);
1317
1318 complex<real_t> val, val_xx, val_xy;
1319 val = 0.25_r * zi * Ho;
1320 val_xx = 0.25_r * zi * (r_xx * Ho_r + r_x * r_x * Ho_rr);
1321 val_xy = 0.25_r * zi * (r_xy * Ho_r + r_x * r_y * Ho_rr);
1322 E[0] = zi / k * (k * k * val + val_xx);
1323 E[1] = zi / k * val_xy;
1324 }
1325 else
1326 {
1327 real_t x0 = X(0) + shift(0);
1328 real_t x1 = X(1) + shift(1);
1329 real_t x2 = X(2) + shift(2);
1330 real_t r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
1331
1332 real_t r_x = x0 / r;
1333 real_t r_y = x1 / r;
1334 real_t r_z = x2 / r;
1335 real_t r_xx = (1.0 / r) * (1.0 - r_x * r_x);
1336 real_t r_yx = -(r_y / r) * r_x;
1337 real_t r_zx = -(r_z / r) * r_x;
1338
1339 complex<real_t> val, val_r, val_rr;
1340 val = exp(zi * k * r) / r;
1341 val_r = val / r * (zi * k * r - 1_r);
1342 val_rr = val / (r * r) * (-k * k * r * r
1343 - 2_r * zi * k * r + 2_r);
1344
1345 complex<real_t> val_xx, val_yx, val_zx;
1346 val_xx = val_rr * r_x * r_x + val_r * r_xx;
1347 val_yx = val_rr * r_x * r_y + val_r * r_yx;
1348 val_zx = val_rr * r_x * r_z + val_r * r_zx;
1349 complex<real_t> alpha = zi * k / 4_r / real_t(M_PI) / k / k;
1350 E[0] = alpha * (k * k * val + val_xx);
1351 E[1] = alpha * val_yx;
1352 E[2] = alpha * val_zx;
1353 }
1354 }
1355 break;
1356
1357 default:
1358 MFEM_ABORT("Should be unreachable");
1359 break;
1360 }
1361}
1362
1364 std::vector<complex<real_t>> &curlE)
1365{
1366 complex<real_t> zi = complex<real_t>(0., 1.);
1367 curlE.resize(dimc);
1368 for (int i = 0; i < dimc; ++i)
1369 {
1370 curlE[i] = 0.0;
1371 }
1372 switch (prob)
1373 {
1374 case plane_wave:
1375 {
1376 std::complex<real_t> pw = exp(zi * omega * (X.Sum()));
1377 if (dim == 3)
1378 {
1379 curlE[0] = 0.0;
1380 curlE[1] = zi * omega * pw;
1381 curlE[2] = -zi * omega * pw;
1382 }
1383 else
1384 {
1385 curlE[0] = -zi * omega * pw;
1386 }
1387 }
1388 break;
1390 {
1391 std::complex<real_t> pw = exp(zi * omega * (X(0)));
1392 curlE[0] = zi * omega * pw;
1393 }
1394 break;
1395 default:
1396 MFEM_ABORT("Should be unreachable");
1397 break;
1398 }
1399}
1400
1402 std::vector<complex<real_t>> &curlcurlE)
1403{
1404 complex<real_t> zi = complex<real_t>(0., 1.);
1405 curlcurlE.resize(dim);
1406 for (int i = 0; i < dim; ++i)
1407 {
1408 curlcurlE[i] = 0.0;
1409 }
1410 switch (prob)
1411 {
1412 case plane_wave:
1413 {
1414 std::complex<real_t> pw = exp(zi * omega * (X.Sum()));
1415 if (dim == 3)
1416 {
1417 curlcurlE[0] = 2_r * omega * omega * pw;
1418 curlcurlE[1] = - omega * omega * pw;
1419 curlcurlE[2] = - omega * omega * pw;
1420 }
1421 else
1422 {
1423 curlcurlE[0] = omega * omega * pw;
1424 curlcurlE[1] = -omega * omega * pw;
1425 }
1426 }
1427 break;
1429 {
1430 std::complex<real_t> pw = exp(zi * omega * (X(0)));
1431 curlcurlE[1] = omega * omega * pw;
1432 }
1433 break;
1434 default:
1435 MFEM_ABORT("Should be unreachable");
1436 break;
1437 }
1438}
1439
1441{
1442 Vector center(dim);
1443 center = 0.5;
1444 real_t r = 0.0;
1445 for (int i = 0; i < dim; ++i)
1446 {
1447 r += pow(x[i] - center[i], 2.);
1448 }
1449 real_t n = 5.0 * omega * sqrt(epsilon * mu) / M_PI;
1450 real_t coeff = pow(n, 2) / M_PI;
1451 real_t alpha = -pow(n, 2) * r;
1452 f = 0.0;
1453 f[0] = -omega * coeff * exp(alpha)/omega;
1454}
Dynamic 2D array using row-major layout.
Definition array.hpp:372
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
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
T & Last()
Return the last element in the array.
Definition array.hpp:802
A class to handle Block diagonal preconditioners in a matrix-free implementation.
Operator & GetDiagonalBlock(int iblock)
Return a reference to block i,i.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, real_t c=1.0)
Add a block op in the block-entry (iblock, jblock).
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
int NumRowBlocks() const
Return the number of row blocks.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
Class for setting up a simple Cartesian PML region.
Definition pml.hpp:19
void SetEpsilonAndMu(real_t epsilon_, real_t mu_)
Definition pml.hpp:65
void SetAttributes(Mesh *mesh_, Array< int > *attrNonPML=nullptr, Array< int > *attrPML=nullptr)
Mark element in the PML region.
Definition pml.cpp:70
void SetOmega(real_t omega_)
Definition pml.hpp:64
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Mimic the action of a complex operator using two real operators.
virtual Operator & imag()
virtual Operator & real()
Real or imaginary part accessor methods.
A coefficient that is constant across space and time.
Integrator for for Nedelec elements.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, const Array< int > &bdr_attr)
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:322
The Auxiliary-space Maxwell Solver in hypre.
Definition hypre.hpp:1845
void SetPrintLevel(int print_lvl)
Definition hypre.cpp:5805
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition hypre.cpp:5111
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Abstract class for hypre's solvers and preconditioners.
Definition hypre.hpp:1162
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
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
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Base class for Matrix Coefficients that optionally depend on time and space.
A matrix coefficient that is constant in space and time.
Matrix coefficient defined as the product of two matrices.
Derived matrix coefficient that has the value of the parent matrix coefficient where it is active and...
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
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10558
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in 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
Arbitrary order H(curl)-trace finite elements defined on the interface between mesh elements (faces,...
Definition fe_coll.hpp:511
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
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 representing the parallel weak formulation. (Convenient for DPG Equations)
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
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
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void MakeRef(FiniteElementSpace *f, real_t *v) override
Make the ParGridFunction reference external data on a new FiniteElementSpace.
void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, const Array< int > &bdr_attr) override
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
Class for parallel meshes.
Definition pmesh.hpp:34
Helper class for ParaView visualization data.
void SetHighOrderOutput(bool high_order_output_)
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
void SetDataFormat(VTKFormat fmt)
Scalar coefficient defined as the product of two scalar coefficients or a scalar and a scalar coeffic...
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:445
Derived coefficient that takes the value of the parent coefficient for the active attributes and is z...
Matrix coefficient defined as a product of a scalar coefficient and a matrix coefficient.
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:922
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
Vector beta
const real_t alpha
Definition ex15.cpp:369
bool exact_known
Definition ex25.cpp:144
prob_type
Definition ex25.cpp:149
int main()
real_t a
Definition lissajous.cpp:41
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
real_t detJ_r_function(const Vector &x, CartesianPML *pml)
PML stretching functions: See https://doi.org/10.1006/jcph.1994.1159.
Definition pml.cpp:160
const int visport
void detJ_Jt_J_inv_r_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:242
float real_t
Definition config.hpp:43
real_t abs_detJ_2_function(const Vector &x, CartesianPML *pml)
Definition pml.cpp:180
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
real_t detJ_i_function(const Vector &x, CartesianPML *pml)
Definition pml.cpp:170
void abs_detJ_Jt_J_inv_2_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:276
void detJ_Jt_J_inv_i_function(const Vector &x, CartesianPML *pml, DenseMatrix &M)
Definition pml.cpp:259
const char vishost[]
void curlE_exact_i(const Vector &x, Vector &curlE_i)
void H_exact_i(const Vector &x, Vector &H_i)
real_t omega
Definition pmaxwell.cpp:189
void curlcurlE_exact_r(const Vector &x, Vector &curlcurlE_r)
int dimc
Definition pmaxwell.cpp:188
void source_function(const Vector &x, Vector &f)
void hatH_exact_r(const Vector &X, Vector &hatH_r)
void hatE_exact_r(const Vector &X, Vector &hatE_r)
void H_exact_r(const Vector &x, Vector &H_r)
real_t hatH_exact_scalar_i(const Vector &X)
void curlH_exact_r(const Vector &x, Vector &curlH_r)
real_t hatH_exact_scalar_r(const Vector &X)
int dim
Definition pmaxwell.cpp:187
real_t mu
Definition pmaxwell.cpp:190
void rhs_func_r(const Vector &x, Vector &J_r)
void hatH_exact_i(const Vector &X, Vector &hatH_i)
void rhs_func_i(const Vector &x, Vector &J_i)
real_t epsilon
Definition pmaxwell.cpp:191
void curlcurlE_exact_i(const Vector &x, Vector &curlcurlE_i)
void E_exact_r(const Vector &x, Vector &E_r)
void maxwell_solution_curlcurl(const Vector &X, std::vector< complex< real_t > > &curlcurlE)
void E_exact_i(const Vector &x, Vector &E_i)
void hatE_exact_i(const Vector &X, Vector &hatE_i)
prob_type prob
Definition pmaxwell.cpp:211
void maxwell_solution_curl(const Vector &X, std::vector< complex< real_t > > &curlE)
prob_type
Definition pmaxwell.cpp:194
@ fichera_oven
Definition pmaxwell.cpp:196
@ plane_wave
Definition pmaxwell.cpp:195
@ pml_plane_wave_scatter
Definition pmaxwell.cpp:198
@ pml_pointsource
Definition pmaxwell.cpp:199
@ pml_general
Definition pmaxwell.cpp:197
void maxwell_solution(const Vector &X, std::vector< complex< real_t > > &E)
void curlE_exact_r(const Vector &x, Vector &curlE_r)
void curlH_exact_i(const Vector &x, Vector &curlH_i)
Helper struct to convert a C++ type to an MPI type.