// Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced // at the Lawrence Livermore National Laboratory. All Rights reserved. See files // LICENSE and NOTICE for details. LLNL-CODE-806117. // // This file is part of the MFEM library. For more information and source code // availability visit https://mfem.org. // // MFEM is free software; you can redistribute it and/or modify it under the // terms of the BSD-3 license. We welcome feedback and contributions, see file // CONTRIBUTING.md for details. // // MFEM Ultraweak DPG Maxwell parallel example // // Compile with: make pmaxwell // // sample run // mpirun -np 4 pmaxwell -m ../../data/star.mesh -o 2 -sref 0 -pref 3 -rnum 0.5 -prob 0 // mpirun -np 4 pmaxwell -m ../../data/inline-quad.mesh -o 3 -sref 0 -pref 3 -rnum 4.8 -sc -prob 0 // mpirun -np 4 pmaxwell -m ../../data/inline-hex.mesh -o 2 -sref 0 -pref 1 -rnum 0.8 -sc -prob 0 // mpirun -np 4 pmaxwell -m ../../data/inline-quad.mesh -o 3 -sref 1 -pref 3 -rnum 4.8 -sc -prob 2 // mpirun -np 4 pmaxwell -o 3 -sref 1 -pref 2 -rnum 11.8 -sc -prob 3 // mpirun -np 4 pmaxwell -o 3 -sref 1 -pref 2 -rnum 9.8 -sc -prob 4 // AMR run. Note that this is a computationally intensive sample run. // We recommend trying it on a large machine with more mpi ranks // mpirun -np 4 pmaxwell -o 3 -sref 0 -pref 15 -prob 1 -theta 0.7 -sc // Description: // This example code demonstrates the use of MFEM to define and solve // the "ultraweak" (UW) DPG formulation for the Maxwell problem // ∇×(1/μ ∇×E) - ω² ϵ E = Ĵ , in Ω // E×n = E₀ , on ∂Ω // It solves the following kinds of problems // 1) Known exact solutions with error convergence rates // a) A manufactured solution problem where E is a plane beam // 2) Fichera "microwave" problem // 3) PML problems // a) Generic PML problem with point source given by the load // b) Plane wave scattering from a square // c) PML problem with a point source prescribed on the boundary // The DPG UW deals with the First Order System // i ω μ H + ∇ × E = 0, in Ω // -i ω ϵ E + ∇ × H = J, in Ω // E × n = E_0, on ∂Ω // Note: Ĵ = -iωJ // The ultraweak-DPG formulation is obtained by integration by parts of both // equations and the introduction of trace unknowns on the mesh skeleton // in 2D // E is vector valued and H is scalar. // (∇ × E, F) = (E, ∇ × F) + < n × E , F> // or (∇ ⋅ AE , F) = (AE, ∇ F) + < AE ⋅ n, F> // where A = [0 1; -1 0]; // E ∈ (L²(Ω))² , H ∈ L²(Ω) // Ê ∈ H^-1/2(Ω)(Γₕ), Ĥ ∈ H^1/2(Γₕ) // i ω μ (H,F) + (E, ∇ × F) + < AÊ, F > = 0, ∀ F ∈ H¹ // -i ω ϵ (E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω) // Ê = E₀ on ∂Ω // ------------------------------------------------------------------------- // | | E | H | Ê | Ĥ | RHS | // ------------------------------------------------------------------------- // | F | (E,∇ × F) | i ω μ (H,F) | < Ê, F > | | | // | | | | | | | // | G | -i ω ϵ (E,G) | (H,∇ × G) | | < Ĥ, G × n > | (J,G) | // where (F,G) ∈ H¹ × H(curl,Ω) // in 3D // E,H ∈ (L^2(Ω))³ // Ê ∈ H_0^1/2(Ω)(curl, Γₕ), Ĥ ∈ H^-1/2(curl, Γₕ) // i ω μ (H,F) + (E,∇ × F) + < Ê, F × n > = 0, ∀ F ∈ H(curl,Ω) // -i ω ϵ (E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω) // Ê × n = E₀ on ∂Ω // ------------------------------------------------------------------------- // | | E | H | Ê | Ĥ | RHS | // ------------------------------------------------------------------------- // | F | (E,∇ × F) | i ω μ (H,F) | < n × Ê, F > | | | // | | | | | | | // | G | -i ω ϵ (E,G) | (H,∇ × G) | | < n × Ĥ, G > | (J,G) | // where (F,G) ∈ H(curl,Ω) × H(curl,Ω) // Here we use the "Adjoint Graph" norm on the test space i.e., // ||(F,G)||²ᵥ = ||A^*(F,G)||² + ||(F,G)||² where A is the // maxwell operator defined by (1) // The PML formulation is // ∇×(1/μ α ∇×E) - ω² ϵ β E = Ĵ , in Ω // E×n = E₀ , on ∂Ω // where α = |J|⁻¹ Jᵀ J (in 2D it's the scalar |J|⁻¹), // β = |J| J⁻¹ J⁻ᵀ, J is the Jacobian of the stretching map // and |J| its determinant. // The first order system reads // i ω μ α⁻¹ H + ∇ × E = 0, in Ω // -i ω ϵ β E + ∇ × H = J, in Ω // E × n = E₀, on ∂Ω // and the ultraweak formulation is // in 2D // E ∈ (L²(Ω))² , H ∈ L²(Ω) // Ê ∈ H^-1/2(Ω)(Γₕ), Ĥ ∈ H^1/2(Γₕ) // i ω μ (α⁻¹ H,F) + (E, ∇ × F) + < AÊ, F > = 0, ∀ F ∈ H¹ // -i ω ϵ (β E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω) // Ê = E₀ on ∂Ω // --------------------------------------------------------------------------------- // | | E | H | Ê | Ĥ | RHS | // --------------------------------------------------------------------------------- // | F | (E,∇ × F) | i ω μ (α⁻¹ H,F) | < Ê, F > | | | // | | | | | | | // | G | -i ω ϵ (β E,G) | (H,∇ × G) | | < Ĥ, G × n > | (J,G) | // where (F,G) ∈ H¹ × H(curl,Ω) // // in 3D // E,H ∈ (L^2(Ω))³ // Ê ∈ H_0^1/2(Ω)(curl, Γ_h), Ĥ ∈ H^-1/2(curl, Γₕ) // i ω μ (α⁻¹ H,F) + (E,∇ × F) + < Ê, F × n > = 0, ∀ F ∈ H(curl,Ω) // -i ω ϵ (β E,G) + (H,∇ × G) + < Ĥ, G × n > = (J,G) ∀ G ∈ H(curl,Ω) // Ê × n = E_0 on ∂Ω // ------------------------------------------------------------------------------- // | | E | H | Ê | Ĥ | RHS | // ------------------------------------------------------------------------------- // | F | ( E,∇ × F) | i ω μ (α⁻¹ H,F) | < n × Ê, F > | | | // | | | | | | | // | G | -iωϵ (β E,G) | (H,∇ × G) | | < n × Ĥ, G > | (J,G) | // where (F,G) ∈ H(curl,Ω) × H(curl,Ω) // For more information see https://doi.org/10.1016/j.camwa.2021.01.017 #include "mfem.hpp" #include "util/pcomplexweakform.hpp" #include "util/pml.hpp" #include "../common/mfem-common.hpp" #include #include using namespace std; using namespace mfem; using namespace mfem::common; void E_exact_r(const Vector &x, Vector & E_r); void E_exact_i(const Vector &x, Vector & E_i); void H_exact_r(const Vector &x, Vector & H_r); void H_exact_i(const Vector &x, Vector & H_i); void rhs_func_r(const Vector &x, Vector & J_r); void rhs_func_i(const Vector &x, Vector & J_i); void curlE_exact_r(const Vector &x, Vector &curlE_r); void curlE_exact_i(const Vector &x, Vector &curlE_i); void curlH_exact_r(const Vector &x,Vector &curlH_r); void curlH_exact_i(const Vector &x,Vector &curlH_i); void curlcurlE_exact_r(const Vector &x, Vector & curlcurlE_r); void curlcurlE_exact_i(const Vector &x, Vector & curlcurlE_i); void hatE_exact_r(const Vector & X, Vector & hatE_r); void hatE_exact_i(const Vector & X, Vector & hatE_i); void hatH_exact_r(const Vector & X, Vector & hatH_r); void hatH_exact_i(const Vector & X, Vector & hatH_i); double hatH_exact_scalar_r(const Vector & X); double hatH_exact_scalar_i(const Vector & X); void maxwell_solution(const Vector & X, std::vector> &E); void maxwell_solution_curl(const Vector & X, std::vector> &curlE); void maxwell_solution_curlcurl(const Vector & X, std::vector> &curlcurlE); void source_function(const Vector &x, Vector & f); int dim; int dimc; double omega; double mu = 1.0; double epsilon = 1.0; enum prob_type { plane_wave, fichera_oven, pml_general, pml_plane_wave_scatter, pml_pointsource }; static const char *enum_str[] = { "plane_wave", "fichera_oven", "pml_general", "pml_plane_wave_scatter", "pml_pointsource" }; prob_type prob; int main(int argc, char *argv[]) { Mpi::Init(); int myid = Mpi::WorldRank(); Hypre::Init(); const char *mesh_file = "../../data/inline-quad.mesh"; int order = 1; int delta_order = 1; double rnum=1.0; double theta = 0.0; bool static_cond = false; int iprob = 0; int sr = 0; int pr = 1; bool exact_known = false; bool with_pml = false; bool visualization = true; bool paraview = false; OptionsParser args(argc, argv); args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use."); args.AddOption(&order, "-o", "--order", "Finite element order (polynomial degree)"); args.AddOption(&rnum, "-rnum", "--number-of-wavelengths", "Number of wavelengths"); args.AddOption(&mu, "-mu", "--permeability", "Permeability of free space (or 1/(spring constant))."); args.AddOption(&epsilon, "-eps", "--permittivity", "Permittivity of free space (or mass constant)."); args.AddOption(&iprob, "-prob", "--problem", "Problem case" " 0: plane wave, 1: Fichera 'oven', " " 2: Generic PML problem with point source given as a load " " 3: Scattering of a plane wave, " " 4: Point source given on the boundary"); args.AddOption(&delta_order, "-do", "--delta-order", "Order enrichment for DPG test space."); args.AddOption(&theta, "-theta", "--theta", "Theta parameter for AMR"); args.AddOption(&sr, "-sref", "--serial-ref", "Number of parallel refinements."); args.AddOption(&pr, "-pref", "--parallel-ref", "Number of parallel refinements."); args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc", "--no-static-condensation", "Enable static condensation."); args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization."); args.AddOption(¶view, "-paraview", "--paraview", "-no-paraview", "--no-paraview", "Enable or disable ParaView visualization."); args.Parse(); if (!args.Good()) { if (myid == 0) { args.PrintUsage(cout); } return 1; } if (iprob > 4) { iprob = 0; } prob = (prob_type)iprob; omega = 2.*M_PI*rnum; if (prob == 0) { exact_known = true; } else if (prob == 1) { mesh_file = "meshes/fichera-waveguide.mesh"; omega = 5.0; rnum = omega/(2.*M_PI); } else if (prob == 2) { with_pml = true; } else { with_pml = true; mesh_file = "meshes/scatter.mesh"; } if (myid == 0) { args.PrintOptions(cout); } Mesh mesh(mesh_file, 1, 1); dim = mesh.Dimension(); MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example"); dimc = (dim == 3) ? 3 : 1; for (int i = 0; i length(dim, 2); length = 0.25; pml = new CartesianPML(&mesh,length); pml->SetOmega(omega); pml->SetEpsilonAndMu(epsilon,mu); } ParMesh pmesh(MPI_COMM_WORLD, mesh); mesh.Clear(); // PML element attribute marker Array attr; Array attrPML; if (pml) { pml->SetAttributes(&pmesh, &attr, &attrPML); } // Define spaces enum TrialSpace { E_space = 0, H_space = 1, hatE_space = 2, hatH_space = 3 }; enum TestSpace { F_space = 0, G_space = 1 }; // L2 space for E FiniteElementCollection *E_fec = new L2_FECollection(order-1,dim); ParFiniteElementSpace *E_fes = new ParFiniteElementSpace(&pmesh,E_fec,dim); // Vector L2 space for H FiniteElementCollection *H_fec = new L2_FECollection(order-1,dim); ParFiniteElementSpace *H_fes = new ParFiniteElementSpace(&pmesh,H_fec, dimc); // H^-1/2 (curl) space for Ê FiniteElementCollection * hatE_fec = nullptr; FiniteElementCollection * hatH_fec = nullptr; FiniteElementCollection * F_fec = nullptr; int test_order = order+delta_order; if (dim == 3) { hatE_fec = new ND_Trace_FECollection(order,dim); hatH_fec = new ND_Trace_FECollection(order,dim); F_fec = new ND_FECollection(test_order, dim); } else { hatE_fec = new RT_Trace_FECollection(order-1,dim); hatH_fec = new H1_Trace_FECollection(order,dim); F_fec = new H1_FECollection(test_order, dim); } ParFiniteElementSpace *hatE_fes = new ParFiniteElementSpace(&pmesh,hatE_fec); ParFiniteElementSpace *hatH_fes = new ParFiniteElementSpace(&pmesh,hatH_fec); FiniteElementCollection * G_fec = new ND_FECollection(test_order, dim); Array trial_fes; Array test_fec; trial_fes.Append(E_fes); trial_fes.Append(H_fes); trial_fes.Append(hatE_fes); trial_fes.Append(hatH_fes); test_fec.Append(F_fec); test_fec.Append(G_fec); // Bilinear form coefficients ConstantCoefficient one(1.0); ConstantCoefficient eps2omeg2(epsilon*epsilon*omega*omega); ConstantCoefficient mu2omeg2(mu*mu*omega*omega); ConstantCoefficient muomeg(mu*omega); ConstantCoefficient negepsomeg(-epsilon*omega); ConstantCoefficient epsomeg(epsilon*omega); ConstantCoefficient negmuomeg(-mu*omega); // for the 2D case DenseMatrix rot_mat(2); rot_mat(0,0) = 0.; rot_mat(0,1) = 1.; rot_mat(1,0) = -1.; rot_mat(1,1) = 0.; MatrixConstantCoefficient rot(rot_mat); ScalarMatrixProductCoefficient epsrot(epsomeg,rot); ScalarMatrixProductCoefficient negepsrot(negepsomeg,rot); Coefficient * epsomeg_cf = nullptr; Coefficient * negepsomeg_cf = nullptr; Coefficient * eps2omeg2_cf = nullptr; Coefficient * muomeg_cf = nullptr; Coefficient * negmuomeg_cf = nullptr; Coefficient * mu2omeg2_cf = nullptr; MatrixCoefficient *epsrot_cf = nullptr; MatrixCoefficient *negepsrot_cf = nullptr; if (pml) { epsomeg_cf = new RestrictedCoefficient(epsomeg,attr); negepsomeg_cf = new RestrictedCoefficient(negepsomeg,attr); eps2omeg2_cf = new RestrictedCoefficient(eps2omeg2,attr); muomeg_cf = new RestrictedCoefficient(muomeg,attr); negmuomeg_cf = new RestrictedCoefficient(negmuomeg,attr); mu2omeg2_cf = new RestrictedCoefficient(mu2omeg2,attr); epsrot_cf = new MatrixRestrictedCoefficient(epsrot,attr); negepsrot_cf = new MatrixRestrictedCoefficient(negepsrot,attr); } else { epsomeg_cf = &epsomeg; negepsomeg_cf = &negepsomeg; eps2omeg2_cf = &eps2omeg2; muomeg_cf = &muomeg; negmuomeg_cf = &negmuomeg; mu2omeg2_cf = &mu2omeg2; epsrot_cf = &epsrot; negepsrot_cf = &negepsrot; } // PML coefficients; PmlCoefficient detJ_r(detJ_r_function,pml); PmlCoefficient detJ_i(detJ_i_function,pml); PmlCoefficient abs_detJ_2(abs_detJ_2_function,pml); PmlMatrixCoefficient detJ_Jt_J_inv_r(dim,detJ_Jt_J_inv_r_function,pml); PmlMatrixCoefficient detJ_Jt_J_inv_i(dim,detJ_Jt_J_inv_i_function,pml); PmlMatrixCoefficient abs_detJ_Jt_J_inv_2(dim,abs_detJ_Jt_J_inv_2_function,pml); ProductCoefficient negmuomeg_detJ_r(negmuomeg,detJ_r); ProductCoefficient negmuomeg_detJ_i(negmuomeg,detJ_i); ProductCoefficient muomeg_detJ_r(muomeg,detJ_r); ProductCoefficient mu2omeg2_detJ_2(mu2omeg2,abs_detJ_2); ScalarMatrixProductCoefficient epsomeg_detJ_Jt_J_inv_i(epsomeg, detJ_Jt_J_inv_i); ScalarMatrixProductCoefficient epsomeg_detJ_Jt_J_inv_r(epsomeg, detJ_Jt_J_inv_r); ScalarMatrixProductCoefficient negepsomeg_detJ_Jt_J_inv_r(negepsomeg, detJ_Jt_J_inv_r); ScalarMatrixProductCoefficient muomeg_detJ_Jt_J_inv_r(muomeg,detJ_Jt_J_inv_r); ScalarMatrixProductCoefficient negmuomeg_detJ_Jt_J_inv_i(negmuomeg, detJ_Jt_J_inv_i); ScalarMatrixProductCoefficient negmuomeg_detJ_Jt_J_inv_r(negmuomeg, detJ_Jt_J_inv_r); ScalarMatrixProductCoefficient mu2omeg2_detJ_Jt_J_inv_2(mu2omeg2, abs_detJ_Jt_J_inv_2); ScalarMatrixProductCoefficient eps2omeg2_detJ_Jt_J_inv_2(eps2omeg2, abs_detJ_Jt_J_inv_2); RestrictedCoefficient negmuomeg_detJ_r_restr(negmuomeg_detJ_r,attrPML); RestrictedCoefficient negmuomeg_detJ_i_restr(negmuomeg_detJ_i,attrPML); RestrictedCoefficient muomeg_detJ_r_restr(muomeg_detJ_r,attrPML); RestrictedCoefficient mu2omeg2_detJ_2_restr(mu2omeg2_detJ_2,attrPML); MatrixRestrictedCoefficient epsomeg_detJ_Jt_J_inv_i_restr( epsomeg_detJ_Jt_J_inv_i,attrPML); MatrixRestrictedCoefficient epsomeg_detJ_Jt_J_inv_r_restr( epsomeg_detJ_Jt_J_inv_r,attrPML); MatrixRestrictedCoefficient negepsomeg_detJ_Jt_J_inv_r_restr( negepsomeg_detJ_Jt_J_inv_r,attrPML); MatrixRestrictedCoefficient muomeg_detJ_Jt_J_inv_r_restr(muomeg_detJ_Jt_J_inv_r, attrPML); MatrixRestrictedCoefficient negmuomeg_detJ_Jt_J_inv_i_restr( negmuomeg_detJ_Jt_J_inv_i,attrPML); MatrixRestrictedCoefficient negmuomeg_detJ_Jt_J_inv_r_restr( negmuomeg_detJ_Jt_J_inv_r,attrPML); MatrixRestrictedCoefficient mu2omeg2_detJ_Jt_J_inv_2_restr( mu2omeg2_detJ_Jt_J_inv_2,attrPML); MatrixRestrictedCoefficient eps2omeg2_detJ_Jt_J_inv_2_restr( eps2omeg2_detJ_Jt_J_inv_2,attrPML); MatrixProductCoefficient * epsomeg_detJ_Jt_J_inv_i_rot = nullptr; MatrixProductCoefficient * epsomeg_detJ_Jt_J_inv_r_rot = nullptr; MatrixProductCoefficient * negepsomeg_detJ_Jt_J_inv_r_rot = nullptr; MatrixRestrictedCoefficient * epsomeg_detJ_Jt_J_inv_i_rot_restr = nullptr; MatrixRestrictedCoefficient * epsomeg_detJ_Jt_J_inv_r_rot_restr = nullptr; MatrixRestrictedCoefficient * negepsomeg_detJ_Jt_J_inv_r_rot_restr = nullptr; if (pml && dim == 2) { epsomeg_detJ_Jt_J_inv_i_rot = new MatrixProductCoefficient( epsomeg_detJ_Jt_J_inv_i, rot); epsomeg_detJ_Jt_J_inv_r_rot = new MatrixProductCoefficient( epsomeg_detJ_Jt_J_inv_r, rot); negepsomeg_detJ_Jt_J_inv_r_rot = new MatrixProductCoefficient( negepsomeg_detJ_Jt_J_inv_r, rot); epsomeg_detJ_Jt_J_inv_i_rot_restr = new MatrixRestrictedCoefficient( *epsomeg_detJ_Jt_J_inv_i_rot, attrPML); epsomeg_detJ_Jt_J_inv_r_rot_restr = new MatrixRestrictedCoefficient( *epsomeg_detJ_Jt_J_inv_r_rot, attrPML); negepsomeg_detJ_Jt_J_inv_r_rot_restr = new MatrixRestrictedCoefficient( *negepsomeg_detJ_Jt_J_inv_r_rot, attrPML); } ParComplexDPGWeakForm * a = new ParComplexDPGWeakForm(trial_fes,test_fec); a->StoreMatrices(); // needed for AMR // (E,∇ × F) a->AddTrialIntegrator(new TransposeIntegrator(new MixedCurlIntegrator(one)), nullptr,TrialSpace::E_space, TestSpace::F_space); // -i ω ϵ (E , G) = i (- ω ϵ E, G) a->AddTrialIntegrator(nullptr, new TransposeIntegrator(new VectorFEMassIntegrator(*negepsomeg_cf)), TrialSpace::E_space,TestSpace::G_space); // (H,∇ × G) a->AddTrialIntegrator(new TransposeIntegrator(new MixedCurlIntegrator(one)), nullptr,TrialSpace::H_space, TestSpace::G_space); // < n×Ĥ ,G> a->AddTrialIntegrator(new TangentTraceIntegrator,nullptr, TrialSpace::hatH_space, TestSpace::G_space); // test integrators // (∇×G ,∇× δG) a->AddTestIntegrator(new CurlCurlIntegrator(one),nullptr, TestSpace::G_space,TestSpace::G_space); // (G,δG) a->AddTestIntegrator(new VectorFEMassIntegrator(one),nullptr, TestSpace::G_space,TestSpace::G_space); if (dim == 3) { // i ω μ (H, F) a->AddTrialIntegrator(nullptr, new TransposeIntegrator( new VectorFEMassIntegrator(*muomeg_cf)), TrialSpace::H_space,TestSpace::F_space); // < n×Ê,F> a->AddTrialIntegrator(new TangentTraceIntegrator,nullptr, TrialSpace::hatE_space, TestSpace::F_space); // test integrators // (∇×F,∇×δF) a->AddTestIntegrator(new CurlCurlIntegrator(one),nullptr, TestSpace::F_space, TestSpace::F_space); // (F,δF) a->AddTestIntegrator(new VectorFEMassIntegrator(one),nullptr, TestSpace::F_space,TestSpace::F_space); // μ^2 ω^2 (F,δF) a->AddTestIntegrator(new VectorFEMassIntegrator(*mu2omeg2_cf),nullptr, TestSpace::F_space, TestSpace::F_space); // -i ω μ (F,∇ × δG) = i (F, -ω μ ∇ × δ G) a->AddTestIntegrator(nullptr,new MixedVectorWeakCurlIntegrator(*negmuomeg_cf), TestSpace::F_space, TestSpace::G_space); // -i ω ϵ (∇ × F, δG) a->AddTestIntegrator(nullptr,new MixedVectorCurlIntegrator(*negepsomeg_cf), TestSpace::F_space, TestSpace::G_space); // i ω μ (∇ × G,δF) a->AddTestIntegrator(nullptr,new MixedVectorCurlIntegrator(*muomeg_cf), TestSpace::G_space, TestSpace::F_space); // i ω ϵ (G, ∇ × δF ) a->AddTestIntegrator(nullptr,new MixedVectorWeakCurlIntegrator(*epsomeg_cf), TestSpace::G_space, TestSpace::F_space); // ϵ^2 ω^2 (G,δG) a->AddTestIntegrator(new VectorFEMassIntegrator(*eps2omeg2_cf),nullptr, TestSpace::G_space, TestSpace::G_space); } else { // i ω μ (H, F) a->AddTrialIntegrator(nullptr,new MixedScalarMassIntegrator(*muomeg_cf), TrialSpace::H_space, TestSpace::F_space); // < n×Ê,F> a->AddTrialIntegrator(new TraceIntegrator,nullptr, TrialSpace::hatE_space, TestSpace::F_space); // test integrators // (∇F,∇δF) a->AddTestIntegrator(new DiffusionIntegrator(one),nullptr, TestSpace::F_space, TestSpace::F_space); // (F,δF) a->AddTestIntegrator(new MassIntegrator(one),nullptr, TestSpace::F_space, TestSpace::F_space); // μ^2 ω^2 (F,δF) a->AddTestIntegrator(new MassIntegrator(*mu2omeg2_cf),nullptr, TestSpace::F_space, TestSpace::F_space); // -i ω μ (F,∇ × δG) = i (F, -ω μ ∇ × δ G) a->AddTestIntegrator(nullptr, new TransposeIntegrator(new MixedCurlIntegrator(*negmuomeg_cf)), TestSpace::F_space, TestSpace::G_space); // -i ω ϵ (∇ × F, δG) = i (- ω ϵ A ∇ F,δG), A = [0 1; -1; 0] a->AddTestIntegrator(nullptr,new MixedVectorGradientIntegrator(*negepsrot_cf), TestSpace::F_space, TestSpace::G_space); // i ω μ (∇ × G,δF) = i (ω μ ∇ × G, δF ) a->AddTestIntegrator(nullptr,new MixedCurlIntegrator(*muomeg_cf), TestSpace::G_space, TestSpace::F_space); // i ω ϵ (G, ∇ × δF ) = i (ω ϵ G, A ∇ δF) = i ( G , ω ϵ A ∇ δF) a->AddTestIntegrator(nullptr, new TransposeIntegrator( new MixedVectorGradientIntegrator(*epsrot_cf)), TestSpace::G_space, TestSpace::F_space); // ϵ^2 ω^2 (G, δG) a->AddTestIntegrator(new VectorFEMassIntegrator(*eps2omeg2_cf),nullptr, TestSpace::G_space, TestSpace::G_space); } if (pml) { //trial integrators // -i ω ϵ (β E , G) = -i ω ϵ ((β_re + i β_im) E, G) // = (ω ϵ β_im E, G) + i (- ω ϵ β_re E, G) a->AddTrialIntegrator( new TransposeIntegrator(new VectorFEMassIntegrator( epsomeg_detJ_Jt_J_inv_i_restr)), new TransposeIntegrator(new VectorFEMassIntegrator( negepsomeg_detJ_Jt_J_inv_r_restr)), TrialSpace::E_space,TestSpace::G_space); if (dim == 3) { //trial integrators // i ω μ (α^-1 H, F) = i ω μ ( (α^-1_re + i α^-1_im) H, F) // = (- ω μ α^-1_im, H,F) + i *(ω μ α^-1_re H, F) a->AddTrialIntegrator( new TransposeIntegrator(new VectorFEMassIntegrator( negmuomeg_detJ_Jt_J_inv_i_restr)), new TransposeIntegrator(new VectorFEMassIntegrator( muomeg_detJ_Jt_J_inv_r_restr)), TrialSpace::H_space, TestSpace::F_space); // test integrators // μ^2 ω^2 (|α|^-2 F,δF) a->AddTestIntegrator( new VectorFEMassIntegrator(mu2omeg2_detJ_Jt_J_inv_2_restr),nullptr, TestSpace::F_space, TestSpace::F_space); // -i ω μ (α^-* F,∇ × δG) = i (F, - ω μ α^-1 ∇ × δ G) // = i (F, - ω μ (α^-1_re + i α^-1_im) ∇ × δ G) // = (F, - ω μ α^-1_im ∇ × δ G) + i (F, - ω μ α^-1_re ∇×δG) a->AddTestIntegrator(new MixedVectorWeakCurlIntegrator( negmuomeg_detJ_Jt_J_inv_i_restr), new MixedVectorWeakCurlIntegrator(negmuomeg_detJ_Jt_J_inv_r_restr), TestSpace::F_space,TestSpace::G_space); // -i ω ϵ (β ∇ × F, δG) = -i ω ϵ ((β_re + i β_im) ∇ × F, δG) // = (ω ϵ β_im ∇ × F, δG) + i (- ω ϵ β_re ∇ × F, δG) a->AddTestIntegrator(new MixedVectorCurlIntegrator( epsomeg_detJ_Jt_J_inv_i_restr), new MixedVectorCurlIntegrator(negepsomeg_detJ_Jt_J_inv_r_restr), TestSpace::F_space,TestSpace::G_space); // i ω μ (α^-1 ∇ × G,δF) = i ω μ ((α^-1_re + i α^-1_im) ∇ × G,δF) // = (- ω μ α^-1_im ∇ × G,δF) + i (ω μ α^-1_re ∇ × G,δF) a->AddTestIntegrator(new MixedVectorCurlIntegrator( negmuomeg_detJ_Jt_J_inv_i_restr), new MixedVectorCurlIntegrator(muomeg_detJ_Jt_J_inv_r_restr), TestSpace::G_space, TestSpace::F_space); // i ω ϵ (β^* G, ∇×δF) = i ω ϵ ( (β_re - i β_im) G, ∇×δF) // = (ω ϵ β_im G, ∇×δF) + i ( ω ϵ β_re G, ∇×δF) a->AddTestIntegrator(new MixedVectorWeakCurlIntegrator( epsomeg_detJ_Jt_J_inv_i_restr), new MixedVectorWeakCurlIntegrator(epsomeg_detJ_Jt_J_inv_r_restr), TestSpace::G_space, TestSpace::F_space); // ϵ^2 ω^2 (|β|^2 G,δG) a->AddTestIntegrator(new VectorFEMassIntegrator( eps2omeg2_detJ_Jt_J_inv_2_restr),nullptr, TestSpace::G_space, TestSpace::G_space); } else { //trial integrators // i ω μ (α^-1 H, F) = i ω μ ( (α^-1_re + i α^-1_im) H, F) // = (- ω μ α^-1_im, H,F) + i *(ω μ α^-1_re H, F) a->AddTrialIntegrator( new MixedScalarMassIntegrator(negmuomeg_detJ_i_restr), new MixedScalarMassIntegrator(muomeg_detJ_r_restr), TrialSpace::H_space, TestSpace::F_space); // test integrators // μ^2 ω^2 (|α|^-2 F,δF) a->AddTestIntegrator(new MassIntegrator(mu2omeg2_detJ_2_restr),nullptr, TestSpace::F_space, TestSpace::F_space); // -i ω μ (α^-* F,∇ × δG) = (F, ω μ α^-1 ∇ × δ G) // =(F, - ω μ α^-1_im ∇ × δ G) + i (F, - ω μ α^-1_re ∇×δG) a->AddTestIntegrator( new TransposeIntegrator(new MixedCurlIntegrator(negmuomeg_detJ_i_restr)), new TransposeIntegrator(new MixedCurlIntegrator(negmuomeg_detJ_r_restr)), TestSpace::F_space, TestSpace::G_space); // -i ω ϵ (β ∇ × F, δG) = i (- ω ϵ β A ∇ F,δG), A = [0 1; -1; 0] // = (ω ϵ β_im A ∇ F, δG) + i (- ω ϵ β_re A ∇ F, δG) a->AddTestIntegrator(new MixedVectorGradientIntegrator( *epsomeg_detJ_Jt_J_inv_i_rot_restr), new MixedVectorGradientIntegrator(*negepsomeg_detJ_Jt_J_inv_r_rot_restr), TestSpace::F_space, TestSpace::G_space); // i ω μ (α^-1 ∇ × G,δF) = i (ω μ α^-1 ∇ × G, δF ) // = (- ω μ α^-1_im ∇ × G,δF) + i (ω μ α^-1_re ∇ × G,δF) a->AddTestIntegrator(new MixedCurlIntegrator(negmuomeg_detJ_i_restr), new MixedCurlIntegrator(muomeg_detJ_r_restr), TestSpace::G_space, TestSpace::F_space); // i ω ϵ (β^* G, ∇ × δF ) = i ( G , ω ϵ β A ∇ δF) // = ( G , ω ϵ β_im A ∇ δF) + i ( G , ω ϵ β_re A ∇ δF) a->AddTestIntegrator( new TransposeIntegrator(new MixedVectorGradientIntegrator( *epsomeg_detJ_Jt_J_inv_i_rot_restr)), new TransposeIntegrator(new MixedVectorGradientIntegrator( *epsomeg_detJ_Jt_J_inv_r_rot_restr)), TestSpace::G_space, TestSpace::F_space); // ϵ^2 ω^2 (|β|^2 G,δG) a->AddTestIntegrator(new VectorFEMassIntegrator( eps2omeg2_detJ_Jt_J_inv_2_restr),nullptr, TestSpace::G_space, TestSpace::G_space); } } // RHS VectorFunctionCoefficient f_rhs_r(dim,rhs_func_r); VectorFunctionCoefficient f_rhs_i(dim,rhs_func_i); VectorFunctionCoefficient f_source(dim,source_function); if (prob == 0) { a->AddDomainLFIntegrator(new VectorFEDomainLFIntegrator(f_rhs_r), new VectorFEDomainLFIntegrator(f_rhs_i), TestSpace::G_space); } else if (prob == 2) { a->AddDomainLFIntegrator(new VectorFEDomainLFIntegrator(f_source),nullptr, TestSpace::G_space); } VectorFunctionCoefficient hatEex_r(dim,hatE_exact_r); VectorFunctionCoefficient hatEex_i(dim,hatE_exact_i); socketstream E_out_r; socketstream H_out_r; if (myid == 0) { std::cout << "\n Ref |" << " Dofs |" << " ω |" ; if (exact_known) { std::cout << " L2 Error |" << " Rate |" ; } std::cout << " Residual |" << " Rate |" << " PCG it |" << endl; std::cout << std::string((exact_known) ? 82 : 60,'-') << endl; } double res0 = 0.; double err0 = 0.; int dof0 = 0; // init to suppress gcc warning Array elements_to_refine; ParGridFunction E_r, E_i, H_r, H_i; ParaViewDataCollection * paraview_dc = nullptr; if (paraview) { paraview_dc = new ParaViewDataCollection(enum_str[prob], &pmesh); paraview_dc->SetPrefixPath("ParaView/Maxwell"); paraview_dc->SetLevelsOfDetail(order); paraview_dc->SetCycle(0); paraview_dc->SetDataFormat(VTKFormat::BINARY); paraview_dc->SetHighOrderOutput(true); paraview_dc->SetTime(0.0); // set the time paraview_dc->RegisterField("E_r",&E_r); paraview_dc->RegisterField("E_i",&E_i); paraview_dc->RegisterField("H_r",&H_r); paraview_dc->RegisterField("H_i",&H_i); } if (static_cond) { a->EnableStaticCondensation(); } for (int it = 0; it<=pr; it++) { a->Assemble(); Array ess_tdof_list; Array ess_bdr; if (pmesh.bdr_attributes.Size()) { ess_bdr.SetSize(pmesh.bdr_attributes.Max()); ess_bdr = 1; hatE_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); if (pml) { ess_bdr = 0; ess_bdr[1] = 1; } } // shift the ess_tdofs for (int j = 0; j < ess_tdof_list.Size(); j++) { ess_tdof_list[j] += E_fes->GetTrueVSize() + H_fes->GetTrueVSize(); } Array offsets(5); offsets[0] = 0; offsets[1] = E_fes->GetVSize(); offsets[2] = H_fes->GetVSize(); offsets[3] = hatE_fes->GetVSize(); offsets[4] = hatH_fes->GetVSize(); offsets.PartialSum(); Vector x(2*offsets.Last()); x = 0.; if (prob != 2) { ParGridFunction hatE_gf_r(hatE_fes, x, offsets[2]); ParGridFunction hatE_gf_i(hatE_fes, x, offsets.Last() + offsets[2]); if (dim == 3) { hatE_gf_r.ProjectBdrCoefficientTangent(hatEex_r, ess_bdr); hatE_gf_i.ProjectBdrCoefficientTangent(hatEex_i, ess_bdr); } else { hatE_gf_r.ProjectBdrCoefficientNormal(hatEex_r, ess_bdr); hatE_gf_i.ProjectBdrCoefficientNormal(hatEex_i, ess_bdr); } } OperatorPtr Ah; Vector X,B; a->FormLinearSystem(ess_tdof_list,x,Ah, X,B); ComplexOperator * Ahc = Ah.As(); BlockOperator * BlockA_r = dynamic_cast(&Ahc->real()); BlockOperator * BlockA_i = dynamic_cast(&Ahc->imag()); int num_blocks = BlockA_r->NumRowBlocks(); Array tdof_offsets(2*num_blocks+1); tdof_offsets[0] = 0; int skip = (static_cond) ? 0 : 2; int k = (static_cond) ? 2 : 0; for (int i=0; iGetTrueVSize(); tdof_offsets[num_blocks+i+1] = trial_fes[i+k]->GetTrueVSize(); } tdof_offsets.PartialSum(); BlockOperator blockA(tdof_offsets); for (int i = 0; iGetBlock(i,j)); blockA.SetBlock(i,j+num_blocks,&BlockA_i->GetBlock(i,j), -1.0); blockA.SetBlock(i+num_blocks,j+num_blocks,&BlockA_r->GetBlock(i,j)); blockA.SetBlock(i+num_blocks,j,&BlockA_i->GetBlock(i,j)); } } X = 0.; BlockDiagonalPreconditioner M(tdof_offsets); if (!static_cond) { HypreBoomerAMG * solver_E = new HypreBoomerAMG((HypreParMatrix &) BlockA_r->GetBlock(0,0)); solver_E->SetPrintLevel(0); solver_E->SetSystemsOptions(dim); HypreBoomerAMG * solver_H = new HypreBoomerAMG((HypreParMatrix &) BlockA_r->GetBlock(1,1)); solver_H->SetPrintLevel(0); solver_H->SetSystemsOptions(dim); M.SetDiagonalBlock(0,solver_E); M.SetDiagonalBlock(1,solver_H); M.SetDiagonalBlock(num_blocks,solver_E); M.SetDiagonalBlock(num_blocks+1,solver_H); } HypreSolver * solver_hatH = nullptr; HypreAMS * solver_hatE = new HypreAMS((HypreParMatrix &)BlockA_r->GetBlock(skip, skip), hatE_fes); solver_hatE->SetPrintLevel(0); if (dim == 2) { solver_hatH = new HypreBoomerAMG((HypreParMatrix &)BlockA_r->GetBlock(skip+1, skip+1)); dynamic_cast(solver_hatH)->SetPrintLevel(0); } else { solver_hatH = new HypreAMS((HypreParMatrix &)BlockA_r->GetBlock(skip+1,skip+1), hatH_fes); dynamic_cast(solver_hatH)->SetPrintLevel(0); } M.SetDiagonalBlock(skip,solver_hatE); M.SetDiagonalBlock(skip+1,solver_hatH); M.SetDiagonalBlock(skip+num_blocks,solver_hatE); M.SetDiagonalBlock(skip+num_blocks+1,solver_hatH); CGSolver cg(MPI_COMM_WORLD); cg.SetRelTol(1e-6); cg.SetMaxIter(10000); cg.SetPrintLevel(0); cg.SetPreconditioner(M); cg.SetOperator(blockA); cg.Mult(B, X); for (int i = 0; iRecoverFEMSolution(X,x); Vector & residuals = a->ComputeResidual(x); double residual = residuals.Norml2(); double maxresidual = residuals.Max(); double globalresidual = residual * residual; MPI_Allreduce(MPI_IN_PLACE,&maxresidual,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD); MPI_Allreduce(MPI_IN_PLACE,&globalresidual,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); globalresidual = sqrt(globalresidual); E_r.MakeRef(E_fes,x, 0); E_i.MakeRef(E_fes,x, offsets.Last()); H_r.MakeRef(H_fes,x, offsets[1]); H_i.MakeRef(H_fes,x, offsets.Last()+offsets[1]); int dofs = 0; for (int i = 0; iGlobalTrueVSize(); } double L2Error = 0.0; double rate_err = 0.0; if (exact_known) { VectorFunctionCoefficient E_ex_r(dim,E_exact_r); VectorFunctionCoefficient E_ex_i(dim,E_exact_i); VectorFunctionCoefficient H_ex_r(dim,H_exact_r); VectorFunctionCoefficient H_ex_i(dim,H_exact_i); double E_err_r = E_r.ComputeL2Error(E_ex_r); double E_err_i = E_i.ComputeL2Error(E_ex_i); double H_err_r = H_r.ComputeL2Error(H_ex_r); double H_err_i = H_i.ComputeL2Error(H_ex_i); L2Error = sqrt( E_err_r*E_err_r + E_err_i*E_err_i + H_err_r*H_err_r + H_err_i*H_err_i ); rate_err = (it) ? dim*log(err0/L2Error)/log((double)dof0/dofs) : 0.0; err0 = L2Error; } double rate_res = (it) ? dim*log(res0/globalresidual)/log(( double)dof0/dofs) : 0.0; res0 = globalresidual; dof0 = dofs; if (myid == 0) { std::ios oldState(nullptr); oldState.copyfmt(std::cout); std::cout << std::right << std::setw(5) << it << " | " << std::setw(10) << dof0 << " | " << std::setprecision(1) << std::fixed << std::setw(4) << 2.0*rnum << " π | " << std::setprecision(3); if (exact_known) { std::cout << std::setw(10) << std::scientific << err0 << " | " << std::setprecision(2) << std::setw(6) << std::fixed << rate_err << " | " ; } std::cout << std::setprecision(3) << std::setw(10) << std::scientific << res0 << " | " << std::setprecision(2) << std::setw(6) << std::fixed << rate_res << " | " << std::setw(6) << std::fixed << num_iter << " | " << std::endl; std::cout.copyfmt(oldState); } if (visualization) { const char * keys = (it == 0 && dim == 2) ? "jRcml\n" : nullptr; char vishost[] = "localhost"; int visport = 19916; VisualizeField(E_out_r,vishost, visport, E_r, "Numerical Electric field (real part)", 0, 0, 500, 500, keys); VisualizeField(H_out_r,vishost, visport, H_r, "Numerical Magnetic field (real part)", 501, 0, 500, 500, keys); } if (paraview) { paraview_dc->SetCycle(it); paraview_dc->SetTime((double)it); paraview_dc->Save(); } if (it == pr) { break; } if (theta > 0.0) { elements_to_refine.SetSize(0); for (int iel = 0; iel theta * maxresidual) { elements_to_refine.Append(iel); } } pmesh.GeneralRefinement(elements_to_refine,1,1); } else { pmesh.UniformRefinement(); } if (pml) { pml->SetAttributes(&pmesh); } for (int i =0; iUpdate(false); } a->Update(); } if (pml && dim == 2) { delete epsomeg_detJ_Jt_J_inv_i_rot; delete epsomeg_detJ_Jt_J_inv_r_rot; delete negepsomeg_detJ_Jt_J_inv_r_rot; delete epsomeg_detJ_Jt_J_inv_i_rot_restr; delete epsomeg_detJ_Jt_J_inv_r_rot_restr; delete negepsomeg_detJ_Jt_J_inv_r_rot_restr; } if (paraview) { delete paraview_dc; } delete a; delete F_fec; delete G_fec; delete hatH_fes; delete hatH_fec; delete hatE_fes; delete hatE_fec; delete H_fec; delete E_fec; delete H_fes; delete E_fes; return 0; } void E_exact_r(const Vector &x, Vector & E_r) { std::vector> E; maxwell_solution(x,E); E_r.SetSize(E.size()); for (unsigned i = 0; i < E.size(); i++) { E_r[i]= E[i].real(); } } void E_exact_i(const Vector &x, Vector & E_i) { std::vector> E; maxwell_solution(x, E); E_i.SetSize(E.size()); for (unsigned i = 0; i < E.size(); i++) { E_i[i]= E[i].imag(); } } void curlE_exact_r(const Vector &x, Vector &curlE_r) { std::vector> curlE; maxwell_solution_curl(x, curlE); curlE_r.SetSize(curlE.size()); for (unsigned i = 0; i < curlE.size(); i++) { curlE_r[i]= curlE[i].real(); } } void curlE_exact_i(const Vector &x, Vector &curlE_i) { std::vector> curlE; maxwell_solution_curl(x, curlE); curlE_i.SetSize(curlE.size()); for (unsigned i = 0; i < curlE.size(); i++) { curlE_i[i]= curlE[i].imag(); } } void curlcurlE_exact_r(const Vector &x, Vector & curlcurlE_r) { std::vector> curlcurlE; maxwell_solution_curlcurl(x, curlcurlE); curlcurlE_r.SetSize(curlcurlE.size()); for (unsigned i = 0; i < curlcurlE.size(); i++) { curlcurlE_r[i]= curlcurlE[i].real(); } } void curlcurlE_exact_i(const Vector &x, Vector & curlcurlE_i) { std::vector> curlcurlE; maxwell_solution_curlcurl(x, curlcurlE); curlcurlE_i.SetSize(curlcurlE.size()); for (unsigned i = 0; i < curlcurlE.size(); i++) { curlcurlE_i[i]= curlcurlE[i].imag(); } } void H_exact_r(const Vector &x, Vector & H_r) { // H = i ∇ × E / ω μ // H_r = - ∇ × E_i / ω μ Vector curlE_i; curlE_exact_i(x,curlE_i); H_r.SetSize(dimc); for (int i = 0; i> &E) { complex zi = complex(0., 1.); E.resize(dim); for (int i = 0; i < dim; ++i) { E[i] = 0.0; } switch (prob) { case plane_wave: { E[0] = exp(zi * omega * (X.Sum())); } break; case pml_plane_wave_scatter: { E[1] = exp(zi * omega * (X(0))); } break; case fichera_oven: { if (abs(X(2) - 3.0) < 1e-10) { E[0] = sin(M_PI*X(1)); } } break; case pml_pointsource: { Vector shift(dim); double k = omega * sqrt(epsilon * mu); shift = -0.5; if (dim == 2) { double x0 = X(0) + shift(0); double x1 = X(1) + shift(1); double r = sqrt(x0 * x0 + x1 * x1); double beta = k * r; // Bessel functions complex Ho, Ho_r, Ho_rr; Ho = jn(0, beta) + zi * yn(0, beta); Ho_r = -k * (jn(1, beta) + zi * yn(1, beta)); Ho_rr = -k * k * (1.0 / beta * (jn(1, beta) + zi * yn(1, beta)) - (jn(2, beta) + zi * yn(2, beta))); // First derivatives double r_x = x0 / r; double r_y = x1 / r; double r_xy = -(r_x / r) * r_y; double r_xx = (1.0 / r) * (1.0 - r_x * r_x); complex val, val_xx, val_xy; val = 0.25 * zi * Ho; val_xx = 0.25 * zi * (r_xx * Ho_r + r_x * r_x * Ho_rr); val_xy = 0.25 * zi * (r_xy * Ho_r + r_x * r_y * Ho_rr); E[0] = zi / k * (k * k * val + val_xx); E[1] = zi / k * val_xy; } else { double x0 = X(0) + shift(0); double x1 = X(1) + shift(1); double x2 = X(2) + shift(2); double r = sqrt(x0 * x0 + x1 * x1 + x2 * x2); double r_x = x0 / r; double r_y = x1 / r; double r_z = x2 / r; double r_xx = (1.0 / r) * (1.0 - r_x * r_x); double r_yx = -(r_y / r) * r_x; double r_zx = -(r_z / r) * r_x; complex val, val_r, val_rr; val = exp(zi * k * r) / r; val_r = val / r * (zi * k * r - 1.0); val_rr = val / (r * r) * (-k * k * r * r - 2.0 * zi * k * r + 2.0); complex val_xx, val_yx, val_zx; val_xx = val_rr * r_x * r_x + val_r * r_xx; val_yx = val_rr * r_x * r_y + val_r * r_yx; val_zx = val_rr * r_x * r_z + val_r * r_zx; complex alpha = zi * k / 4.0 / M_PI / k / k; E[0] = alpha * (k * k * val + val_xx); E[1] = alpha * val_yx; E[2] = alpha * val_zx; } } break; default: MFEM_ABORT("Should be unreachable"); break; } } void maxwell_solution_curl(const Vector & X, std::vector> &curlE) { complex zi = complex(0., 1.); curlE.resize(dimc); for (int i = 0; i < dimc; ++i) { curlE[i] = 0.0; } switch (prob) { case plane_wave: { std::complex pw = exp(zi * omega * (X.Sum())); if (dim == 3) { curlE[0] = 0.0; curlE[1] = zi * omega * pw; curlE[2] = -zi * omega * pw; } else { curlE[0] = -zi * omega * pw; } } break; case pml_plane_wave_scatter: { std::complex pw = exp(zi * omega * (X(0))); curlE[0] = zi * omega * pw; } break; default: MFEM_ABORT("Should be unreachable"); break; } } void maxwell_solution_curlcurl(const Vector & X, std::vector> &curlcurlE) { complex zi = complex(0., 1.); curlcurlE.resize(dim); for (int i = 0; i < dim; ++i) { curlcurlE[i] = 0.0;; } switch (prob) { case plane_wave: { std::complex pw = exp(zi * omega * (X.Sum())); if (dim == 3) { curlcurlE[0] = 2.0 * omega * omega * pw; curlcurlE[1] = - omega * omega * pw; curlcurlE[2] = - omega * omega * pw; } else { curlcurlE[0] = omega * omega * pw; curlcurlE[1] = -omega * omega * pw; } } break; case pml_plane_wave_scatter: { std::complex pw = exp(zi * omega * (X(0))); curlcurlE[1] = omega * omega * pw; } break; default: MFEM_ABORT("Should be unreachable"); break; } } void source_function(const Vector &x, Vector &f) { Vector center(dim); center = 0.5; double r = 0.0; for (int i = 0; i < dim; ++i) { r += pow(x[i] - center[i], 2.); } double n = 5.0 * omega * sqrt(epsilon * mu) / M_PI; double coeff = pow(n, 2) / M_PI; double alpha = -pow(n, 2) * r; f = 0.0; f[0] = -omega * coeff * exp(alpha)/omega; }