MFEM  v4.6.0 Finite element discretization library
maxwell.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
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 example for Maxwell
13 //
14 // Compile with: make maxwell
15 //
16 // Sample runs
17 // maxwell -m ../../data/inline-tri.mesh -ref 4 -o 1 -rnum 1.0
18 // maxwell -m ../../data/amr-quad.mesh -ref 3 -o 2 -rnum 1.6 -sc
19 // maxwell -m ../../data/inline-quad.mesh -ref 2 -o 3 -rnum 4.2 -sc
20 // maxwell -m ../../data/inline-hex.mesh -ref 1 -o 2 -sc -rnum 1.0
21
22 // Description:
23 // This example code demonstrates the use of MFEM to define and solve
24 // the "ultraweak" (UW) DPG formulation for the (indefinite) Maxwell problem
25
26 // ∇×(1/μ ∇×E) - ω² ϵ E = Ĵ , in Ω
27 // E×n = E₀, on ∂Ω
28
29 // It solves a problem with a manufactured solution E_exact being a plane wave
30 // in the x-component and zero in y (and z) component.
31 // This example computes and prints out convergence rates for the L² error.
32
33 // The DPG UW deals with the First Order System
34 // i ω μ H + ∇ × E = 0, in Ω
35 // -i ω ϵ E + ∇ × H = J, in Ω (1)
36 // E × n = E₀, on ∂Ω
37
38 // Note: Ĵ = -iωJ
39 // in 2D
40 // E is vector valued and H is scalar.
41 // (∇ × E, F) = (E, ∇ × F) + < n × E , F>
42 // or (∇ ⋅ AE , F) = (AE, ∇ F) + < AE ⋅ n, F>
43 // where A = [0 1; -1 0];
44
45 // E ∈ (L²(Ω))² , H ∈ L²(Ω)
46 // Ê ∈ H^-1/2(Ω)(Γₕ), Ĥ ∈ H^1/2(Γₕ)
47 // i ω μ (H,F) + (E, ∇ × F) + < AÊ, F > = 0, ∀ F ∈ H¹
48 // -i ω ϵ (E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω)
49 // Ê = E₀ on ∂Ω
50 // -------------------------------------------------------------------------
51 // | | E | H | Ê | Ĥ | RHS |
52 // -------------------------------------------------------------------------
53 // | F | (E,∇ × F) | i ω μ (H,F) | < Ê, F > | | |
54 // | | | | | | |
55 // | G | -i ω ϵ (E,G) | (H,∇ × G) | | < Ĥ, G × n > | (J,G) |
56 // where (F,G) ∈ H¹ × H(curl,Ω)
57
58 // in 3D
59 // E,H ∈ ((L²(Ω)))³
60 // Ê ∈ H_0^1/2(Ω)(curl, Γₕ), Ĥ ∈ H^-1/2(curl, Γₕ)
61 // i ω μ (H,F) + (E,∇ × F) + < Ê, F × n > = 0, ∀ F ∈ H(curl,Ω)
62 // -i ω ϵ (E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω)
63 // Ê × n = E_0 on ∂Ω
64 // -------------------------------------------------------------------------
65 // | | E | H | Ê | Ĥ | RHS |
66 // -------------------------------------------------------------------------
67 // | F | (E,∇ × F) | i ω μ (H,F) | < n × Ê, F > | | |
68 // | | | | | | |
69 // | G | -i ω ϵ (E,G) | (H,∇ × G) | | < n × Ĥ, G > | (J,G) |
70 // where (F,G) ∈ H(curl,Ω) × H(curl,Ω)
71
72 // Here we use the "Adjoint Graph" norm on the test space i.e.,
73 // ||(F,G)||²ᵥ = ||A^*(F,G)||² + ||(F,G)||² where A is the
74 // Maxwell operator defined by (1)
75
77
78 #include "mfem.hpp"
79 #include "util/complexweakform.hpp"
80 #include "../common/mfem-common.hpp"
81 #include <fstream>
82 #include <iostream>
83
84 using namespace std;
85 using namespace mfem;
86 using namespace mfem::common;
87
88 void maxwell_solution(const Vector & X, std::vector<complex<double>> &E);
89
90 void maxwell_solution_curl(const Vector & X,
91  std::vector<complex<double>> &curlE);
92
93 void maxwell_solution_curlcurl(const Vector & X,
94  std::vector<complex<double>> &curlcurlE);
95
96 void E_exact_r(const Vector &x, Vector & E_r);
97 void E_exact_i(const Vector &x, Vector & E_i);
98
99 void H_exact_r(const Vector &x, Vector & H_r);
100 void H_exact_i(const Vector &x, Vector & H_i);
101
102 void curlE_exact_r(const Vector &x, Vector &curlE_r);
103 void curlE_exact_i(const Vector &x, Vector &curlE_i);
104 void curlH_exact_r(const Vector &x,Vector &curlH_r);
105 void curlH_exact_i(const Vector &x,Vector &curlH_i);
106
107 void curlcurlE_exact_r(const Vector &x, Vector & curlcurlE_r);
108 void curlcurlE_exact_i(const Vector &x, Vector & curlcurlE_i);
109
110 void hatE_exact_r(const Vector & X, Vector & hatE_r);
111 void hatE_exact_i(const Vector & X, Vector & hatE_i);
112
113 void hatH_exact_r(const Vector & X, Vector & hatH_r);
114 void hatH_exact_i(const Vector & X, Vector & hatH_i);
115
116 double hatH_exact_scalar_r(const Vector & X);
117 double hatH_exact_scalar_i(const Vector & X);
118
119 void rhs_func_r(const Vector &x, Vector & J_r);
120 void rhs_func_i(const Vector &x, Vector & J_i);
121
122 int dim;
123 int dimc;
124 double omega;
125 double mu = 1.0;
126 double epsilon = 1.0;
127
128 int main(int argc, char *argv[])
129 {
130  const char *mesh_file = "../../data/inline-quad.mesh";
131  int order = 1;
132  int delta_order = 1;
133  bool visualization = true;
134  double rnum=1.0;
135  int ref = 0;
136  bool static_cond = false;
137
138  OptionsParser args(argc, argv);
140  "Mesh file to use.");
142  "Finite element order (polynomial degree)");
144  "--no-visualization",
145  "Enable or disable GLVis visualization.");
147  "Number of wavelengths");
149  "Permeability of free space (or 1/(spring constant)).");
151  "Permittivity of free space (or mass constant).");
153  "Order enrichment for DPG test space.");
155  "Number of serial refinements.");
157  "--no-static-condensation", "Enable static condensation.");
158  args.Parse();
159  if (!args.Good())
160  {
161  args.PrintUsage(cout);
162  return 1;
163  }
164  args.PrintOptions(cout);
165
166  omega = 2.*M_PI*rnum;
167
168  Mesh mesh(mesh_file, 1, 1);
169  dim = mesh.Dimension();
170  MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
171
172  dimc = (dim == 3) ? 3 : 1;
173  int test_order = order+delta_order;
174
175  // Define spaces
176  enum TrialSpace
177  {
178  E_space = 0,
179  H_space = 1,
180  hatE_space = 2,
181  hatH_space = 3
182  };
183  enum TestSpace
184  {
185  F_space = 0,
186  G_space = 1
187  };
188  // L2 space for E
189  FiniteElementCollection *E_fec = new L2_FECollection(order-1,dim);
190  FiniteElementSpace *E_fes = new FiniteElementSpace(&mesh,E_fec,dim);
191
192  // Vector L2 space for H
193  FiniteElementCollection *H_fec = new L2_FECollection(order-1,dim);
194  FiniteElementSpace *H_fes = new FiniteElementSpace(&mesh,H_fec, dimc);
195
196  // H^-1/2 (curl) space for Ê
197  FiniteElementCollection * hatE_fec = nullptr;
198  FiniteElementCollection * hatH_fec = nullptr;
199  FiniteElementCollection * F_fec = nullptr;
200  if (dim == 3)
201  {
202  hatE_fec = new ND_Trace_FECollection(order,dim);
203  hatH_fec = new ND_Trace_FECollection(order,dim);
204  F_fec = new ND_FECollection(test_order, dim);
205  }
206  else
207  {
208  hatE_fec = new RT_Trace_FECollection(order-1,dim);
209  hatH_fec = new H1_Trace_FECollection(order,dim);
210  F_fec = new H1_FECollection(test_order, dim);
211  }
212  FiniteElementSpace *hatE_fes = new FiniteElementSpace(&mesh,hatE_fec);
213  FiniteElementSpace *hatH_fes = new FiniteElementSpace(&mesh,hatH_fec);
214  FiniteElementCollection * G_fec = new ND_FECollection(test_order, dim);
215
216  // Coefficients
217  ConstantCoefficient one(1.0);
219  ConstantCoefficient mu2omeg2(mu*mu*omega*omega);
220  ConstantCoefficient muomeg(mu*omega);
221  ConstantCoefficient negepsomeg(-epsilon*omega);
223  ConstantCoefficient negmuomeg(-mu*omega);
224
225  DenseMatrix rot_mat(2);
226  rot_mat(0,0) = 0.; rot_mat(0,1) = 1.;
227  rot_mat(1,0) = -1.; rot_mat(1,1) = 0.;
228  MatrixConstantCoefficient rot(rot_mat);
229  ScalarMatrixProductCoefficient epsrot(epsomeg,rot);
230  ScalarMatrixProductCoefficient negepsrot(negepsomeg,rot);
231
234
235  trial_fes.Append(E_fes);
236  trial_fes.Append(H_fes);
237  trial_fes.Append(hatE_fes);
238  trial_fes.Append(hatH_fes);
239
240  test_fec.Append(F_fec);
241  test_fec.Append(G_fec);
242
243  ComplexDPGWeakForm * a = new ComplexDPGWeakForm(trial_fes,test_fec);
244  a->StoreMatrices();
245
246  // (E,∇ × F)
248  nullptr,TrialSpace::E_space,TestSpace::F_space);
249  // -i ω ϵ (E , G)
251  new TransposeIntegrator(new VectorFEMassIntegrator(negepsomeg)),
252  TrialSpace::E_space,TestSpace::G_space);
253  // (H,∇ × G)
255  nullptr, TrialSpace::H_space,TestSpace::G_space);
256
257  if (dim == 3)
258  {
259  // i ω μ (H, F)
262  TrialSpace::H_space,TestSpace::F_space);
263  // < n×Ê,F>
265  TrialSpace::hatE_space,TestSpace::F_space);
266  }
267  else
268  {
269  // i ω μ (H, F)
271  TrialSpace::H_space,TestSpace::F_space);
272  // <Ê,F>
274  TrialSpace::hatE_space, TestSpace::F_space);
275  }
276  // < n×Ĥ ,G>
278  TrialSpace::hatH_space, TestSpace::G_space);
279
280  // test integrators for the adjoint graph norm on the test space
281  // (∇×G ,∇× δG)
283  TestSpace::G_space,TestSpace::G_space);
284  // (G,δG)
286  TestSpace::G_space, TestSpace::G_space);
287
288  if (dim == 3)
289  {
290  // (∇×F,∇×δF)
292  TestSpace::F_space, TestSpace::F_space);
293  // (F,δF)
295  TestSpace::F_space,TestSpace::F_space);
296  // μ^2 ω^2 (F,δF)
298  TestSpace::F_space, TestSpace::F_space);
299  // -i ω μ (F,∇ × δG) = (F, ω μ ∇ × δ G)
301  TestSpace::F_space, TestSpace::G_space);
302  // -i ω ϵ (∇ × F, δG)
304  TestSpace::F_space, TestSpace::G_space);
305  // i ω μ (∇ × G,δF)
307  TestSpace::G_space, TestSpace::F_space);
308  // i ω ϵ (G, ∇ × δF )
310  TestSpace::G_space, TestSpace::F_space);
311  // ϵ^2 ω^2 (G,δG)
313  TestSpace::G_space, TestSpace::G_space);
314  }
315  else
316  {
317  // (∇F,∇δF)
319  TestSpace::F_space, TestSpace::F_space);
320  // (F,δF)
322  TestSpace::F_space, TestSpace::F_space);
323  // μ^2 ω^2 (F,δF)
325  TestSpace::F_space, TestSpace::F_space);
326  // -i ω μ (F,∇ × δG) = i (F, -ω μ ∇ × δ G)
328  new TransposeIntegrator(new MixedCurlIntegrator(negmuomeg)),
329  TestSpace::F_space, TestSpace::G_space);
330  // -i ω ϵ (∇ × F, δG) = i (- ω ϵ A ∇ F,δG), A = [0 1; -1; 0]
332  TestSpace::F_space, TestSpace::G_space);
333  // i ω μ (∇ × G,δF) = i (ω μ ∇ × G, δF )
335  TestSpace::G_space, TestSpace::F_space);
336  // i ω ϵ (G, ∇ × δF ) = i (ω ϵ G, A ∇ δF) = i ( G , ω ϵ A ∇ δF)
339  TestSpace::G_space, TestSpace::F_space);
340  // ϵ^2 ω^2 (G,δG)
342  TestSpace::G_space, TestSpace::G_space);
343  }
344
345  // RHS
349  new VectorFEDomainLFIntegrator(f_rhs_i),
350  TestSpace::G_space);
351
354
355  socketstream E_out_r;
356  socketstream E_out_i;
357
358  double err0 = 0.;
359  int dof0 = 0; // init to suppress gcc warning
360
361  std::cout << "\n Ref |"
362  << " Dofs |"
363  << " ω |"
364  << " L2 Error |"
365  << " Rate |"
366  << " PCG it |" << endl;
367  std::cout << std::string(60,'-')
368  << endl;
369
370  for (int it = 0; it<=ref; it++)
371  {
372  if (static_cond) { a->EnableStaticCondensation(); }
373  a->Assemble();
374
375  Array<int> ess_tdof_list;
376  Array<int> ess_bdr;
377  if (mesh.bdr_attributes.Size())
378  {
379  ess_bdr.SetSize(mesh.bdr_attributes.Max());
380  ess_bdr = 1;
381  hatE_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
382  }
383
384  // shift the ess_tdofs
385  for (int j = 0; j < ess_tdof_list.Size(); j++)
386  {
387  ess_tdof_list[j] += E_fes->GetTrueVSize() + H_fes->GetTrueVSize();
388  }
389
390  Array<int> offsets(5);
391  offsets[0] = 0;
392  offsets[1] = E_fes->GetVSize();
393  offsets[2] = H_fes->GetVSize();
394  offsets[3] = hatE_fes->GetVSize();
395  offsets[4] = hatH_fes->GetVSize();
396  offsets.PartialSum();
397
398  Vector x(2*offsets.Last());
399  x = 0.;
400
401  GridFunction hatE_gf_r(hatE_fes, x, offsets[2]);
402  GridFunction hatE_gf_i(hatE_fes, x, offsets.Last() + offsets[2]);
403  if (dim == 3)
404  {
405  hatE_gf_r.ProjectBdrCoefficientTangent(hatEex_r, ess_bdr);
406  hatE_gf_i.ProjectBdrCoefficientTangent(hatEex_i, ess_bdr);
407  }
408  else
409  {
410  hatE_gf_r.ProjectBdrCoefficientNormal(hatEex_r, ess_bdr);
411  hatE_gf_i.ProjectBdrCoefficientNormal(hatEex_i, ess_bdr);
412  }
413  OperatorPtr Ah;
414  Vector X,B;
415  a->FormLinearSystem(ess_tdof_list,x,Ah, X,B);
416
417  ComplexOperator * Ahc = Ah.As<ComplexOperator>();
418
419  BlockMatrix * A_r = dynamic_cast<BlockMatrix *>(&Ahc->real());
420  BlockMatrix * A_i = dynamic_cast<BlockMatrix *>(&Ahc->imag());
421
422  int num_blocks = A_r->NumRowBlocks();
423  Array<int> tdof_offsets(2*num_blocks+1);
424  tdof_offsets[0] = 0;
425  int k = (static_cond) ? 2 : 0;
426  for (int i=0; i<num_blocks; i++)
427  {
428  tdof_offsets[i+1] = trial_fes[i+k]->GetTrueVSize();
429  tdof_offsets[num_blocks+i+1] = trial_fes[i+k]->GetTrueVSize();
430  }
431  tdof_offsets.PartialSum();
432
433  BlockOperator A(tdof_offsets);
434  for (int i = 0; i<num_blocks; i++)
435  {
436  for (int j = 0; j<num_blocks; j++)
437  {
438  A.SetBlock(i,j,&A_r->GetBlock(i,j));
439  A.SetBlock(i,j+num_blocks,&A_i->GetBlock(i,j), -1.0);
440  A.SetBlock(i+num_blocks,j+num_blocks,&A_r->GetBlock(i,j));
441  A.SetBlock(i+num_blocks,j,&A_i->GetBlock(i,j));
442  }
443  }
444
445  BlockDiagonalPreconditioner M(tdof_offsets);
446  M.owns_blocks = 1;
447  for (int i = 0; i<num_blocks; i++)
448  {
449  M.SetDiagonalBlock(i, new GSSmoother((SparseMatrix&)A_r->GetBlock(i,i)));
450  M.SetDiagonalBlock(num_blocks+i, new GSSmoother((SparseMatrix&)A_r->GetBlock(i,
451  i)));
452  }
453
454  CGSolver cg;
455  cg.SetRelTol(1e-10);
456  cg.SetMaxIter(2000);
457  cg.SetPrintLevel(0);
458  cg.SetPreconditioner(M);
459  cg.SetOperator(A);
460  cg.Mult(B, X);
461
462  a->RecoverFEMSolution(X,x);
463
464  GridFunction E_r(E_fes, x, 0);
465  GridFunction E_i(E_fes, x, offsets.Last());
466
469
470  GridFunction H_r(H_fes, x, offsets[1]);
471  GridFunction H_i(H_fes, x, offsets.Last() + offsets[1]);
472
475
476  int dofs = 0;
477  for (int i = 0; i<trial_fes.Size(); i++)
478  {
479  dofs += trial_fes[i]->GetTrueVSize();
480  }
481
482  double E_err_r = E_r.ComputeL2Error(E_ex_r);
483  double E_err_i = E_i.ComputeL2Error(E_ex_i);
484  double H_err_r = H_r.ComputeL2Error(H_ex_r);
485  double H_err_i = H_i.ComputeL2Error(H_ex_i);
486
487  double L2Error = sqrt(E_err_r*E_err_r + E_err_i*E_err_i
488  + H_err_r*H_err_r + H_err_i*H_err_i);
489
490  double rate_err = (it) ? dim*log(err0/L2Error)/log((double)dof0/dofs) : 0.0;
491
492  err0 = L2Error;
493  dof0 = dofs;
494
495  std::ios oldState(nullptr);
496  oldState.copyfmt(std::cout);
497  std::cout << std::right << std::setw(5) << it << " | "
498  << std::setw(10) << dof0 << " | "
499  << std::setprecision(1) << std::fixed
500  << std::setw(4) << 2*rnum << " π | "
501  << std::setprecision(3)
502  << std::setw(10) << std::scientific << err0 << " | "
503  << std::setprecision(2)
504  << std::setw(6) << std::fixed << rate_err << " | "
505  << std::setw(6) << std::fixed << cg.GetNumIterations() << " | "
506  << std::endl;
507  std::cout.copyfmt(oldState);
508
509  if (visualization)
510  {
511  const char * keys = (it == 0 && dim == 2) ? "jRcml\n" : nullptr;
512  char vishost[] = "localhost";
513  int visport = 19916;
514  VisualizeField(E_out_r,vishost, visport, E_r,
515  "Numerical Electric field (real part)", 0, 0, 500, 500, keys);
516  VisualizeField(E_out_i,vishost, visport, E_i,
517  "Numerical Electric field (imaginary part)", 501, 0, 500, 500, keys);
518  }
519
520  if (it == ref)
521  {
522  break;
523  }
524
525  mesh.UniformRefinement();
526  for (int i =0; i<trial_fes.Size(); i++)
527  {
528  trial_fes[i]->Update(false);
529  }
530  a->Update();
531  }
532
533  delete a;
534  delete F_fec;
535  delete G_fec;
536  delete hatH_fes;
537  delete hatH_fec;
538  delete hatE_fes;
539  delete hatE_fec;
540  delete H_fec;
541  delete E_fec;
542  delete H_fes;
543  delete E_fes;
544
545  return 0;
546 }
547
548 void E_exact_r(const Vector &x, Vector & E_r)
549 {
550  std::vector<std::complex<double>> E;
551  maxwell_solution(x, E);
552  E_r.SetSize(E.size());
553  for (unsigned i = 0; i < E.size(); i++)
554  {
555  E_r[i]= E[i].real();
556  }
557 }
558
559 void E_exact_i(const Vector &x, Vector & E_i)
560 {
561  std::vector<std::complex<double>> E;
562  maxwell_solution(x, E);
563  E_i.SetSize(E.size());
564  for (unsigned i = 0; i < E.size(); i++)
565  {
566  E_i[i]= E[i].imag();
567  }
568 }
569
570 void curlE_exact_r(const Vector &x, Vector &curlE_r)
571 {
572  std::vector<std::complex<double>> curlE;
573  maxwell_solution_curl(x, curlE);
574  curlE_r.SetSize(curlE.size());
575  for (unsigned i = 0; i < curlE.size(); i++)
576  {
577  curlE_r[i]= curlE[i].real();
578  }
579 }
580
581 void curlE_exact_i(const Vector &x, Vector &curlE_i)
582 {
583  std::vector<std::complex<double>> curlE;
584  maxwell_solution_curl(x, curlE);
585  curlE_i.SetSize(curlE.size());
586  for (unsigned i = 0; i < curlE.size(); i++)
587  {
588  curlE_i[i]= curlE[i].imag();
589  }
590 }
591
592 void curlcurlE_exact_r(const Vector &x, Vector & curlcurlE_r)
593 {
594  std::vector<std::complex<double>> curlcurlE;
595  maxwell_solution_curlcurl(x, curlcurlE);
596  curlcurlE_r.SetSize(curlcurlE.size());
597  for (unsigned i = 0; i < curlcurlE.size(); i++)
598  {
599  curlcurlE_r[i]= curlcurlE[i].real();
600  }
601 }
602
603 void curlcurlE_exact_i(const Vector &x, Vector & curlcurlE_i)
604 {
605  std::vector<std::complex<double>> curlcurlE;
606  maxwell_solution_curlcurl(x, curlcurlE);
607  curlcurlE_i.SetSize(curlcurlE.size());
608  for (unsigned i = 0; i < curlcurlE.size(); i++)
609  {
610  curlcurlE_i[i]= curlcurlE[i].imag();
611  }
612 }
613
614 void H_exact_r(const Vector &x, Vector & H_r)
615 {
616  // H = i ∇ × E / ω μ
617  // H_r = - ∇ × E_i / ω μ
618  Vector curlE_i;
619  curlE_exact_i(x,curlE_i);
620  H_r.SetSize(dimc);
621  for (int i = 0; i<dimc; i++)
622  {
623  H_r(i) = - curlE_i(i) / (omega * mu);
624  }
625 }
626
627 void H_exact_i(const Vector &x, Vector & H_i)
628 {
629  // H = i ∇ × E / ω μ
630  // H_i = ∇ × E_r / ω μ
631  Vector curlE_r;
632  curlE_exact_r(x,curlE_r);
633  H_i.SetSize(dimc);
634  for (int i = 0; i<dimc; i++)
635  {
636  H_i(i) = curlE_r(i) / (omega * mu);
637  }
638 }
639
640 void curlH_exact_r(const Vector &x,Vector &curlH_r)
641 {
642  // ∇ × H_r = - ∇ × ∇ × E_i / ω μ
643  Vector curlcurlE_i;
644  curlcurlE_exact_i(x,curlcurlE_i);
645  curlH_r.SetSize(dim);
646  for (int i = 0; i<dim; i++)
647  {
648  curlH_r(i) = -curlcurlE_i(i) / (omega * mu);
649  }
650 }
651
652 void curlH_exact_i(const Vector &x,Vector &curlH_i)
653 {
654  // ∇ × H_i = ∇ × ∇ × E_r / ω μ
655  Vector curlcurlE_r;
656  curlcurlE_exact_r(x,curlcurlE_r);
657  curlH_i.SetSize(dim);
658  for (int i = 0; i<dim; i++)
659  {
660  curlH_i(i) = curlcurlE_r(i) / (omega * mu);
661  }
662 }
663
664 void hatE_exact_r(const Vector & x, Vector & hatE_r)
665 {
666  if (dim == 3)
667  {
668  E_exact_r(x,hatE_r);
669  }
670  else
671  {
672  Vector E_r;
673  E_exact_r(x,E_r);
674  hatE_r.SetSize(hatE_r.Size());
675  // rotate E_hat
676  hatE_r[0] = E_r[1];
677  hatE_r[1] = -E_r[0];
678  }
679 }
680
681 void hatE_exact_i(const Vector & x, Vector & hatE_i)
682 {
683  if (dim == 3)
684  {
685  E_exact_i(x,hatE_i);
686  }
687  else
688  {
689  Vector E_i;
690  E_exact_i(x,E_i);
691  hatE_i.SetSize(hatE_i.Size());
692  // rotate E_hat
693  hatE_i[0] = E_i[1];
694  hatE_i[1] = -E_i[0];
695  }
696 }
697
698 void hatH_exact_r(const Vector & x, Vector & hatH_r)
699 {
700  H_exact_r(x,hatH_r);
701 }
702
703 void hatH_exact_i(const Vector & x, Vector & hatH_i)
704 {
705  H_exact_i(x,hatH_i);
706 }
707
708 double hatH_exact_scalar_r(const Vector & x)
709 {
710  Vector hatH_r;
711  H_exact_r(x,hatH_r);
712  return hatH_r[0];
713 }
714
715 double hatH_exact_scalar_i(const Vector & x)
716 {
717  Vector hatH_i;
718  H_exact_i(x,hatH_i);
719  return hatH_i[0];
720 }
721
722 // J = -i ω ϵ E + ∇ × H
723 // J_r + iJ_i = -i ω ϵ (E_r + i E_i) + ∇ × (H_r + i H_i)
724 void rhs_func_r(const Vector &x, Vector & J_r)
725 {
726  // J_r = ω ϵ E_i + ∇ × H_r
727  Vector E_i, curlH_r;
728  E_exact_i(x,E_i);
729  curlH_exact_r(x,curlH_r);
730  J_r.SetSize(dim);
731  for (int i = 0; i<dim; i++)
732  {
733  J_r(i) = omega * epsilon * E_i(i) + curlH_r(i);
734  }
735 }
736
737 void rhs_func_i(const Vector &x, Vector & J_i)
738 {
739  // J_i = - ω ϵ E_r + ∇ × H_i
740  Vector E_r, curlH_i;
741  E_exact_r(x,E_r);
742  curlH_exact_i(x,curlH_i);
743  J_i.SetSize(dim);
744  for (int i = 0; i<dim; i++)
745  {
746  J_i(i) = -omega * epsilon * E_r(i) + curlH_i(i);
747  }
748 }
749
750 void maxwell_solution(const Vector & X, std::vector<complex<double>> &E)
751 {
752  E.resize(dim);
753  std::complex<double> zi(0,1);
754  std::complex<double> pw = exp(-zi * omega * (X.Sum()));
755  E[0] = pw;
756  E[1] = 0.0;
757  if (dim == 3) { E[2] = 0.0; }
758 }
759
761  std::vector<complex<double>> &curlE)
762 {
763  curlE.resize(dimc);
764  std::complex<double> zi(0,1);
765  std::complex<double> pw = exp(-zi * omega * (X.Sum()));
766  if (dim == 3)
767  {
768  curlE[0] = 0.0;
769  curlE[1] = -zi * omega * pw;
770  curlE[2] = zi * omega * pw;
771  }
772  else
773  {
774  curlE[0] = zi * omega * pw;
775  }
776 }
777
779  std::vector<complex<double>> &curlcurlE)
780 {
781  curlcurlE.resize(dim);
782  std::complex<double> zi(0,1);
783  std::complex<double> pw = exp(-zi * omega * (X.Sum()));
784  if (dim == 3)
785  {
786  curlcurlE[0] = 2.0 * omega * omega * pw;
787  curlcurlE[1] = - omega * omega * pw;
788  curlcurlE[2] = - omega * omega * pw;
789  }
790  else
791  {
792  curlcurlE[0] = omega * omega * pw;
793  curlcurlE[1] = - omega * omega * pw ;
794  }
795 }
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:2786
A matrix coefficient that is constant in space and time.
double mu
Definition: maxwell.cpp:125
int visport
void maxwell_solution_curlcurl(const Vector &X, std::vector< complex< double >> &curlcurlE)
Definition: maxwell.cpp:778
Definition: solvers.hpp:493
virtual Operator & real()
Real or imaginary part accessor methods.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
Mimic the action of a complex operator using two real operators.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Integrator for (curl u, curl v) for Nedelec elements.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void curlH_exact_i(const Vector &x, Vector &curlH_i)
Definition: maxwell.cpp:652
void curlE_exact_r(const Vector &x, Vector &curlE_r)
Definition: maxwell.cpp:570
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:718
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void curlH_exact_r(const Vector &x, Vector &curlH_r)
Definition: maxwell.cpp:640
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:587
void hatH_exact_i(const Vector &X, Vector &hatH_i)
Definition: maxwell.cpp:703
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Matrix coefficient defined as a product of a scalar coefficient and a matrix coefficient.
STL namespace.
void curlE_exact_i(const Vector &x, Vector &curlE_i)
Definition: maxwell.cpp:581
double hatH_exact_scalar_i(const Vector &X)
Definition: maxwell.cpp:715
void hatE_exact_r(const Vector &X, Vector &hatE_r)
Definition: maxwell.cpp:664
Data type for Gauss-Seidel smoother of sparse matrix.
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Definition: gridfunc.cpp:2694
int NumRowBlocks() const
Return the number of row blocks.
Definition: blockmatrix.hpp:46
void H_exact_r(const Vector &x, Vector &H_r)
Definition: maxwell.cpp:614
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
Class representing the DPG weak formulation for complex valued systems (see the class DPGWeakForm)...
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
char vishost[]
Data type sparse matrix.
Definition: sparsemat.hpp:50
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
void hatE_exact_i(const Vector &X, Vector &hatE_i)
Definition: maxwell.cpp:681
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition: solvers.hpp:250
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition: fe_coll.hpp:319
Arbitrary order H(curl)-trace finite elements defined on the interface between mesh elements (faces...
Definition: fe_coll.hpp:498
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:1286
int main(int argc, char *argv[])
Definition: maxwell.cpp:127
void maxwell_solution_curl(const Vector &X, std::vector< complex< double >> &curlE)
Definition: maxwell.cpp:760
void curlcurlE_exact_i(const Vector &x, Vector &curlcurlE_i)
Definition: maxwell.cpp:603
A general vector function coefficient.
void hatH_exact_r(const Vector &X, Vector &hatH_r)
Definition: maxwell.cpp:698
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
void E_exact_i(const Vector &x, Vector &E_i)
Definition: maxwell.cpp:559
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
virtual Operator & imag()
void E_exact_r(const Vector &x, Vector &E_r)
Definition: maxwell.cpp:548
int dimc
Definition: maxwell.cpp:123
double a
Definition: lissajous.cpp:41
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
void H_exact_i(const Vector &x, Vector &H_i)
Definition: maxwell.cpp:627
double epsilon
Definition: maxwell.cpp:126
void curlcurlE_exact_r(const Vector &x, Vector &curlcurlE_r)
Definition: maxwell.cpp:592
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
Definition: gridfunc.cpp:2769
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
Definition: blockmatrix.cpp:82
int dim
Definition: maxwell.cpp:122
void rhs_func_i(const Vector &x, Vector &J_i)
Definition: maxwell.cpp:737
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
T & Last()
Return the last element in the array.
Definition: array.hpp:792
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:346
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:454
Vector data type.
Definition: vector.hpp:58
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition: fe_coll.hpp:434
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
double omega
Definition: maxwell.cpp:124
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
double hatH_exact_scalar_r(const Vector &X)
Definition: maxwell.cpp:708
A class to handle Block systems in a matrix-free implementation.
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
void maxwell_solution(const Vector &X, std::vector< complex< double >> &E)
Definition: maxwell.cpp:750
void rhs_func_r(const Vector &x, Vector &J_r)
Definition: maxwell.cpp:724