step-81.h

Go to the documentation of this file.

,
 *                              types::material_id material)
 *     {
 *       return (material == 1 ? epsilon_1 : epsilon_2);
 *     }
 *   
 *     template <int dim>
 *     std::complex<double> Parameters<dim>::mu_inv(const Point<dim> & /*x*/,
 *                                                  types::material_id material)
 *     {
 *       return (material == 1 ? mu_inv_1 : mu_inv_2);
 *     }
 *   
 *     template <int dim>
 *     typename Parameters<dim>::rank2_type
 *     Parameters<dim>::sigma(const Point<dim> & /*x*/,
 *                            types::material_id left,
 *                            types::material_id right)
 *     {
 *       return (left == right ? rank2_type() : sigma_tensor);
 *     }
 *   
 *     template <int dim>
 *     typename Parameters<dim>::rank1_type
 *     Parameters<dim>::J_a(const Point<dim> &point, types::material_id /*id*/)
 *     {
 *       rank1_type J_a;
 *       const auto distance = (dipole_position - point).norm() / dipole_radius;
 *       if (distance > 1.)
 *         return J_a;
 *       double scale = std::cos(distance * numbers::PI / 2.) *
 *                      std::cos(distance * numbers::PI / 2.) /
 *                      (numbers::PI / 2. - 2. / numbers::PI) / dipole_radius /
 *                      dipole_radius;
 *       J_a = dipole_strength * dipole_orientation * scale;
 *       return J_a;
 *     }
 *   
 * @endcode
 * 
 * 
 * <a name="step_81-PerfectlyMatchedLayerClass"></a> 
 * <h4>PerfectlyMatchedLayer Class</h4>
 * The PerfectlyMatchedLayer class inherits ParameterAcceptor as well. It
 * implements the transformation matrices used to modify the permittivity
 * and permeability tensors supplied from the Parameters class. The
 * actual transformation of the material tensors will be done in the
 * assembly loop. The radii and the strength of the PML is specified, and
 * the coefficients will be modified using transformation matrices within
 * the PML region. The radii and strength of the PML are editable through
 * a .prm file. The rotation function @f$T_{exer}@f$ is the same as
 * introduced in the perfectly matched layer section of the introduction.
 * Similarly, the matrices A, B and C are defined as follows
 * @f[
 * A = T_{e_xe_r}^{-1}
 * \text{diag}\left(\frac{1}{\bar{d}^2},\frac{1}{d\bar{d}}\right)T_{e_xe_r},\qquad
 * B = T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{d}\right)T_{e_xe_r},\qquad
 * C = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{d}},\frac{1}{d}\right)
 * T_{e_xe_r}.\qquad
 * @f]
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     class PerfectlyMatchedLayer : public ParameterAcceptor
 *     {
 *     public:
 *       static_assert(dim == 2,
 *                     "The perfectly matched layer is only implemented in 2d.");
 *   
 *       Parameters<dim> parameters;
 *   
 *       using rank1_type = Tensor<1, dim, std::complex<double>>;
 *   
 *       using rank2_type = Tensor<2, dim, std::complex<double>>;
 *   
 *       PerfectlyMatchedLayer();
 *   
 *       std::complex<double> d(const Point<dim> point);
 *   
 *       std::complex<double> d_bar(const Point<dim> point);
 *   
 *   
 *       rank2_type rotation(std::complex<double> d_1,
 *                           std::complex<double> d_2,
 *                           Point<dim>           point);
 *   
 *       rank2_type a_matrix(const Point<dim> point);
 *   
 *       rank2_type b_matrix(const Point<dim> point);
 *   
 *       rank2_type c_matrix(const Point<dim> point);
 *   
 *     private:
 *       double inner_radius;
 *       double outer_radius;
 *       double strength;
 *     };
 *   
 *   
 *     template <int dim>
 *     PerfectlyMatchedLayer<dim>::PerfectlyMatchedLayer()
 *       : ParameterAcceptor("PerfectlyMatchedLayer")
 *     {
 *       inner_radius = 12.;
 *       add_parameter("inner radius",
 *                     inner_radius,
 *                     "inner radius of the PML shell");
 *       outer_radius = 20.;
 *       add_parameter("outer radius",
 *                     outer_radius,
 *                     "outer radius of the PML shell");
 *       strength = 8.;
 *       add_parameter("strength", strength, "strength of the PML");
 *     }
 *   
 *   
 *     template <int dim>
 *     typename std::complex<double>
 *     PerfectlyMatchedLayer<dim>::d(const Point<dim> point)
 *     {
 *       const auto radius = point.norm();
 *       if (radius > inner_radius)
 *         {
 *           const double s =
 *             strength * ((radius - inner_radius) * (radius - inner_radius)) /
 *             ((outer_radius - inner_radius) * (outer_radius - inner_radius));
 *           return {1.0, s};
 *         }
 *       else
 *         {
 *           return 1.0;
 *         }
 *     }
 *   
 *   
 *     template <int dim>
 *     typename std::complex<double>
 *     PerfectlyMatchedLayer<dim>::d_bar(const Point<dim> point)
 *     {
 *       const auto radius = point.norm();
 *       if (radius > inner_radius)
 *         {
 *           const double s_bar =
 *             strength / 3. *
 *             ((radius - inner_radius) * (radius - inner_radius) *
 *              (radius - inner_radius)) /
 *             (radius * (outer_radius - inner_radius) *
 *              (outer_radius - inner_radius));
 *           return {1.0, s_bar};
 *         }
 *       else
 *         {
 *           return 1.0;
 *         }
 *     }
 *   
 *   
 *     template <int dim>
 *     typename PerfectlyMatchedLayer<dim>::rank2_type
 *     PerfectlyMatchedLayer<dim>::rotation(std::complex<double> d_1,
 *                                          std::complex<double> d_2,
 *                                          Point<dim>           point)
 *     {
 *       rank2_type result;
 *       result[0][0] = point[0] * point[0] * d_1 + point[1] * point[1] * d_2;
 *       result[0][1] = point[0] * point[1] * (d_1 - d_2);
 *       result[1][0] = point[0] * point[1] * (d_1 - d_2);
 *       result[1][1] = point[1] * point[1] * d_1 + point[0] * point[0] * d_2;
 *       return result;
 *     }
 *   
 *   
 *     template <int dim>
 *     typename PerfectlyMatchedLayer<dim>::rank2_type
 *     PerfectlyMatchedLayer<dim>::a_matrix(const Point<dim> point)
 *     {
 *       const auto d     = this->d(point);
 *       const auto d_bar = this->d_bar(point);
 *       return invert(rotation(d * d, d * d_bar, point)) *
 *              rotation(d * d, d * d_bar, point);
 *     }
 *   
 *   
 *     template <int dim>
 *     typename PerfectlyMatchedLayer<dim>::rank2_type
 *     PerfectlyMatchedLayer<dim>::b_matrix(const Point<dim> point)
 *     {
 *       const auto d     = this->d(point);
 *       const auto d_bar = this->d_bar(point);
 *       return invert(rotation(d, d_bar, point)) * rotation(d, d_bar, point);
 *     }
 *   
 *   
 *     template <int dim>
 *     typename PerfectlyMatchedLayer<dim>::rank2_type
 *     PerfectlyMatchedLayer<dim>::c_matrix(const Point<dim> point)
 *     {
 *       const auto d     = this->d(point);
 *       const auto d_bar = this->d_bar(point);
 *       return invert(rotation(1. / d_bar, 1. / d, point)) *
 *              rotation(1. / d_bar, 1. / d, point);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_81-MaxwellClass"></a> 
 * <h4>Maxwell Class</h4>
 * At this point we are ready to declare all the major building blocks of
 * the finite element program which consists of the usual setup and
 * assembly routines. Most of the structure has already been introduced
 * in previous tutorial programs. The Maxwell class also holds private
 * instances of the Parameters and PerfectlyMatchedLayers classes
 * introduced above. The default values of these parameters are set to
 * show us a standing wave with absorbing boundary conditions and a PML.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     class Maxwell : public ParameterAcceptor
 *     {
 *     public:
 *       Maxwell();
 *       void run();
 *   
 *     private:
 *       /* run time parameters */
 *       double       scaling;
 *       unsigned int refinements;
 *       unsigned int fe_order;
 *       unsigned int quadrature_order;
 *       bool         absorbing_boundary;
 *   
 *       void parse_parameters_callback();
 *       void make_grid();
 *       void setup_system();
 *       void assemble_system();
 *       void solve();
 *       void output_results();
 *   
 *       Parameters<dim>            parameters;
 *       PerfectlyMatchedLayer<dim> perfectly_matched_layer;
 *   
 *       Triangulation<dim> triangulation;
 *       DoFHandler<dim>    dof_handler;
 *   
 *       std::unique_ptr<FiniteElement<dim>> fe;
 *   
 *       AffineConstraints<double> constraints;
 *       SparsityPattern           sparsity_pattern;
 *       SparseMatrix<double>      system_matrix;
 *       Vector<double>            solution;
 *       Vector<double>            system_rhs;
 *     };
 *   
 * @endcode
 * 
 * 
 * <a name="step_81-ClassTemplateDefinitionsandImplementation"></a> 
 * <h3>Class Template Definitions and Implementation</h3>
 *   
 
 * 
 * 
 * <a name="step_81-TheConstructor"></a> 
 * <h4>The Constructor</h4>
 * The Constructor simply consists of default initialization a number of
 * discretization parameters (such as the domain size, mesh refinement,
 * and the order of finite elements and quadrature) and declaring a
 * corresponding entry via ParameterAcceptor::add_parameter(). All of
 * these can be modified by editing the .prm file. Absorbing boundary
 * conditions can be controlled with the absorbing_boundary boolean. If
 * absorbing boundary conditions are disabled we simply enforce
 * homogeneous Dirichlet conditions on the tangential component of the
 * electric field. In the context of time-harmonic Maxwell's equations
 * these are also known as perfectly conducting boundary conditions.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     Maxwell<dim>::Maxwell()
 *       : ParameterAcceptor("Maxwell")
 *       , dof_handler(triangulation)
 *     {
 *       ParameterAcceptor::parse_parameters_call_back.connect(
 *         [&]() { parse_parameters_callback(); });
 *   
 *       scaling = 20;
 *       add_parameter("scaling", scaling, "scale of the hypercube geometry");
 *   
 *       refinements = 8;
 *       add_parameter("refinements",
 *                     refinements,
 *                     "number of refinements of the geometry");
 *   
 *       fe_order = 0;
 *       add_parameter("fe order", fe_order, "order of the finite element space");
 *   
 *       quadrature_order = 1;
 *       add_parameter("quadrature order",
 *                     quadrature_order,
 *                     "order of the quadrature");
 *   
 *       absorbing_boundary = true;
 *       add_parameter("absorbing boundary condition",
 *                     absorbing_boundary,
 *                     "use absorbing boundary conditions?");
 *     }
 *   
 *   
 *     template <int dim>
 *     void Maxwell<dim>::parse_parameters_callback()
 *     {
 *       fe = std::make_unique<FESystem<dim>>(FE_NedelecSZ<dim>(fe_order), 2);
 *     }
 *   
 * @endcode
 * 
 * The Maxwell::make_grid() routine creates the mesh for the
 * computational domain which in our case is a scaled square domain.
 * Additionally, a material interface is introduced by setting the
 * material id of the upper half (@f$y>0@f$) to 1 and of the lower half
 * (@f$y<0@f$) of the computational domain to 2.
 * We are using a block decomposition into real and imaginary matrices
 * for the solution matrices. More details on this are available
 * under the Results section.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void Maxwell<dim>::make_grid()
 *     {
 *       GridGenerator::hyper_cube(triangulation, -scaling, scaling);
 *       triangulation.refine_global(refinements);
 *   
 *       if (absorbing_boundary)
 *         {
 *           for (auto &face : triangulation.active_face_iterators())
 *             if (face->at_boundary())
 *               face->set_boundary_id(1);
 *         };
 *   
 *       for (auto &cell : triangulation.active_cell_iterators())
 *         if (cell->center()[1] > 0.)
 *           cell->set_material_id(1);
 *         else
 *           cell->set_material_id(2);
 *   
 *   
 *       std::cout << "Number of active cells: " << triangulation.n_active_cells()
 *                 << std::endl;
 *     }
 *   
 * @endcode
 * 
 * The Maxwell::setup_system() routine follows the usual routine of
 * enumerating all the degrees of freedom and setting up the matrix and
 * vector objects to hold the system data. Enumerating is done by using
 * DoFHandler::distribute_dofs().
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void Maxwell<dim>::setup_system()
 *     {
 *       dof_handler.distribute_dofs(*fe);
 *       std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << std::endl;
 *   
 *       solution.reinit(dof_handler.n_dofs());
 *       system_rhs.reinit(dof_handler.n_dofs());
 *   
 *       constraints.clear();
 *   
 *       DoFTools::make_hanging_node_constraints(dof_handler, constraints);
 *   
 *       VectorTools::project_boundary_values_curl_conforming_l2(
 *         dof_handler,
 *         0, /* real part */
 *         Functions::ZeroFunction<dim>(2 * dim),
 *         0, /* boundary id */
 *         constraints);
 *       VectorTools::project_boundary_values_curl_conforming_l2(
 *         dof_handler,
 *         dim, /* imaginary part */
 *         Functions::ZeroFunction<dim>(2 * dim),
 *         0, /* boundary id */
 *         constraints);
 *   
 *       constraints.close();
 *   
 *       DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
 *       DoFTools::make_sparsity_pattern(dof_handler,
 *                                       dsp,
 *                                       constraints,
 *                                       /* keep_constrained_dofs = */ true);
 *       sparsity_pattern.copy_from(dsp);
 *       system_matrix.reinit(sparsity_pattern);
 *     }
 *   
 * @endcode
 * 
 * This is a helper function that takes the tangential component of a tensor.
 * 
 * @code
 *     template <int dim>
 *     DEAL_II_ALWAYS_INLINE inline Tensor<1, dim, std::complex<double>>
 *     tangential_part(const Tensor<1, dim, std::complex<double>> &tensor,
 *                     const Tensor<1, dim>                       &normal)
 *     {
 *       auto result = tensor;
 *       result[0]   = normal[1] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
 *       result[1]   = -normal[0] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
 *       return result;
 *     }
 *   
 *   
 * @endcode
 * 
 * Assemble the stiffness matrix and the right-hand side:
 * \f{align*}{
 * A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_j) \cdot
 * (\nabla\times\bar{\varphi}_i)\text{d}x
 * - \int_\Omega \varepsilon_r\varphi_j \cdot \bar{\varphi}_i\text{d}x
 * - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_j)_T) \cdot
 * (\bar{\varphi}_i)_T\text{do}x
 * - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_j)_T) \cdot
 * (\nabla\times(\bar{\varphi}_i)_T)\text{d}x, \f} \f{align}{
 * F_i = i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x - \int_\Omega
 * \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x.
 * \f}
 * In addition, we will be modifying the coefficients if the position of the
 * cell is within the PML region.
 * 
 
 * 
 * 
 * @code
 *     template <int dim>
 *     void Maxwell<dim>::assemble_system()
 *     {
 *       const QGauss<dim>     quadrature_formula(quadrature_order);
 *       const QGauss<dim - 1> face_quadrature_formula(quadrature_order);
 *   
 *       FEValues<dim, dim>     fe_values(*fe,
 *                                    quadrature_formula,
 *                                    update_values | update_gradients |
 *                                      update_quadrature_points |
 *                                      update_JxW_values);
 *       FEFaceValues<dim, dim> fe_face_values(*fe,
 *                                             face_quadrature_formula,
 *                                             update_values | update_gradients |
 *                                               update_quadrature_points |
 *                                               update_normal_vectors |
 *                                               update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe->dofs_per_cell;
 *   
 *       const unsigned int n_q_points      = quadrature_formula.size();
 *       const unsigned int n_face_q_points = face_quadrature_formula.size();
 *   
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *       Vector<double>     cell_rhs(dofs_per_cell);
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 * @endcode
 * 
 * Next, let us assemble on the interior of the domain on the left hand
 * side. So we are computing
 * \f{align*}{
 * \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot
 * (\nabla\times\bar{\varphi}_j)\text{d}x
 * -
 * \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x
 * \f}
 * and
 * \f{align}{
 * i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x
 * - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i})
 * \text{d}x.
 * \f}
 * In doing so, we need test functions @f$\varphi_i@f$ and @f$\varphi_j@f$, and the
 * curl of these test variables. We must be careful with the signs of the
 * imaginary parts of these complex test variables. Moreover, we have a
 * conditional that changes the parameters if the cell is in the PML region.
 * 
 * @code
 *       const FEValuesExtractors::Vector real_part(0);
 *       const FEValuesExtractors::Vector imag_part(dim);
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *   
 *           cell_matrix = 0.;
 *           cell_rhs    = 0.;
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *           const auto id = cell->material_id();
 *   
 *           const auto &quadrature_points = fe_values.get_quadrature_points();
 *   
 *           for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
 *             {
 *               const Point<dim> &position = quadrature_points[q_point];
 *   
 *               auto       mu_inv  = parameters.mu_inv(position, id);
 *               auto       epsilon = parameters.epsilon(position, id);
 *               const auto J_a     = parameters.J_a(position, id);
 *   
 *               const auto A = perfectly_matched_layer.a_matrix(position);
 *               const auto B = perfectly_matched_layer.b_matrix(position);
 *               const auto d = perfectly_matched_layer.d(position);
 *   
 *               mu_inv  = mu_inv / d;
 *               epsilon = invert(A) * epsilon * invert(B);
 *   
 *               for (const auto i : fe_values.dof_indices())
 *                 {
 *                   constexpr std::complex<double> imag{0., 1.};
 *   
 *                   const auto phi_i =
 *                     fe_values[real_part].value(i, q_point) -
 *                     imag * fe_values[imag_part].value(i, q_point);
 *                   const auto curl_phi_i =
 *                     fe_values[real_part].curl(i, q_point) -
 *                     imag * fe_values[imag_part].curl(i, q_point);
 *   
 *                   const auto rhs_value =
 *                     (imag * (J_a * phi_i)) * fe_values.JxW(q_point);
 *                   cell_rhs(i) += rhs_value.real();
 *   
 *                   for (const auto j : fe_values.dof_indices())
 *                     {
 *                       const auto phi_j =
 *                         fe_values[real_part].value(j, q_point) +
 *                         imag * fe_values[imag_part].value(j, q_point);
 *                       const auto curl_phi_j =
 *                         fe_values[real_part].curl(j, q_point) +
 *                         imag * fe_values[imag_part].curl(j, q_point);
 *   
 *                       const auto temp = (((mu_inv * curl_phi_j) * curl_phi_i) -
 *                                          ((epsilon * phi_j) * phi_i)) *
 *                                         fe_values.JxW(q_point);
 *                       cell_matrix(i, j) += temp.real();
 *                     }
 *                 }
 *             }
 *   
 * @endcode
 * 
 * Now we assemble the face and the boundary. The following loops will
 * assemble
 * \f{align*}{
 * - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot
 * (\bar{\varphi}_j)_T\text{do}x \f} and \f{align}{
 * - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T)
 * \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x,
 * \f}
 * respectively. The test variables and the PML are implemented
 * similarly as the domain.
 *         
 
 * 
 * If we are at the domain boundary @f$\partial\Omega@f$ and absorbing
 * boundary conditions are set (<code>id == 1</code>) we assemble
 * the corresponding boundary term:
 *         
 
 * 
 * 
 * @code
 *           const FEValuesExtractors::Vector real_part(0);
 *           const FEValuesExtractors::Vector imag_part(dim);
 *           for (const auto &face : cell->face_iterators())
 *             {
 *               if (face->at_boundary())
 *                 {
 *                   const auto id = face->boundary_id();
 *                   if (id == 1)
 *                     {
 *                       fe_face_values.reinit(cell, face);
 *   
 *                       for (unsigned int q_point = 0; q_point < n_face_q_points;
 *                            ++q_point)
 *                         {
 *                           const auto &position = quadrature_points[q_point];
 *   
 *                           auto mu_inv  = parameters.mu_inv(position, id);
 *                           auto epsilon = parameters.epsilon(position, id);
 *   
 *                           const auto A =
 *                             perfectly_matched_layer.a_matrix(position);
 *                           const auto B =
 *                             perfectly_matched_layer.b_matrix(position);
 *                           const auto d = perfectly_matched_layer.d(position);
 *   
 *                           mu_inv  = mu_inv / d;
 *                           epsilon = invert(A) * epsilon * invert(B);
 *   
 *                           const auto normal =
 *                             fe_face_values.normal_vector(q_point);
 *   
 *                           for (const auto i : fe_face_values.dof_indices())
 *                             {
 *                               constexpr std::complex<double> imag{0., 1.};
 *   
 *                               const auto phi_i =
 *                                 fe_face_values[real_part].value(i, q_point) -
 *                                 imag *
 *                                   fe_face_values[imag_part].value(i, q_point);
 *                               const auto phi_i_T = tangential_part(phi_i, normal);
 *   
 *                               for (const auto j : fe_face_values.dof_indices())
 *                                 {
 *                                   const auto phi_j =
 *                                     fe_face_values[real_part].value(j, q_point) +
 *                                     imag *
 *                                       fe_face_values[imag_part].value(j, q_point);
 *                                   const auto phi_j_T =
 *                                     tangential_part(phi_j, normal) *
 *                                     fe_face_values.JxW(q_point);
 *   
 *                                   const auto prod      = mu_inv * epsilon;
 *                                   const auto sqrt_prod = prod;
 *   
 *                                   const auto temp =
 *                                     -imag * ((sqrt_prod * phi_j_T) * phi_i_T);
 *                                   cell_matrix(i, j) += temp.real();
 *                                 } /* j */
 *                             }     /* i */
 *                         }         /* q_point */
 *                     }
 *                 }
 *               else
 *                 {
 * @endcode
 * 
 * We are on an interior face:
 * 
 * @code
 *                   const auto face_index = cell->face_iterator_to_index(face);
 *   
 *                   const auto id1 = cell->material_id();
 *                   const auto id2 = cell->neighbor(face_index)->material_id();
 *   
 *                   if (id1 == id2)
 *                     continue; /* skip this face */
 *   
 *                   fe_face_values.reinit(cell, face);
 *   
 *                   for (unsigned int q_point = 0; q_point < n_face_q_points;
 *                        ++q_point)
 *                     {
 *                       const auto &position = quadrature_points[q_point];
 *   
 *                       auto sigma = parameters.sigma(position, id1, id2);
 *   
 *                       const auto B = perfectly_matched_layer.b_matrix(position);
 *                       const auto C = perfectly_matched_layer.c_matrix(position);
 *                       sigma        = invert(C) * sigma * invert(B);
 *   
 *                       const auto normal = fe_face_values.normal_vector(q_point);
 *   
 *                       for (const auto i : fe_face_values.dof_indices())
 *                         {
 *                           constexpr std::complex<double> imag{0., 1.};
 *   
 *                           const auto phi_i =
 *                             fe_face_values[real_part].value(i, q_point) -
 *                             imag * fe_face_values[imag_part].value(i, q_point);
 *                           const auto phi_i_T = tangential_part(phi_i, normal);
 *   
 *                           for (const auto j : fe_face_values.dof_indices())
 *                             {
 *                               const auto phi_j =
 *                                 fe_face_values[real_part].value(j, q_point) +
 *                                 imag *
 *                                   fe_face_values[imag_part].value(j, q_point);
 *                               const auto phi_j_T = tangential_part(phi_j, normal);
 *   
 *                               const auto temp = -imag *
 *                                                 ((sigma * phi_j_T) * phi_i_T) *
 *                                                 fe_face_values.JxW(q_point);
 *                               cell_matrix(i, j) += temp.real();
 *                             } /* j */
 *                         }     /* i */
 *                     }         /* q_point */
 *                 }
 *             }
 *   
 *           constraints.distribute_local_to_global(
 *             cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
 *         }
 *     }
 *   
 * @endcode
 * 
 * We use a direct solver from the SparseDirectUMFPACK to solve the system
 * 
 * @code
 *     template <int dim>
 *     void Maxwell<dim>::solve()
 *     {
 *       SparseDirectUMFPACK A_direct;
 *       A_direct.initialize(system_matrix);
 *       A_direct.vmult(solution, system_rhs);
 *     }
 *   
 * @endcode
 * 
 * The output is written into a vtk file with 4 components
 * 
 * @code
 *     template <int dim>
 *     void Maxwell<dim>::output_results()
 *     {
 *       DataOut<2> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution,
 *                                {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"});
 *       data_out.build_patches();
 *       const std::string filename = "solution.vtk";
 *       std::ofstream     output(filename);
 *       data_out.write_vtk(output);
 *       std::cout << "Output written to " << filename << std::endl;
 *     }
 *   
 *   
 *     template <int dim>
 *     void Maxwell<dim>::run()
 *     {
 *       make_grid();
 *       setup_system();
 *       assemble_system();
 *       solve();
 *       output_results();
 *     }
 *   
 *   } // namespace Step81
 *   
 * @endcode
 * 
 * The following main function calls the class @ref step_81 "step-81"(), initializes the
 * ParameterAcceptor, and calls the run() function.
 * 
 
 * 
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         using namespace dealii;
 *   
 *         Step81::Maxwell<2> maxwell_2d;
 *         ParameterAcceptor::initialize("parameters.prm");
 *         maxwell_2d.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *     return 0;
 *   }
 * @endcode
<a name="step_81-Results"></a><h1>Results</h1>
 
 
The solution is written to a .vtk file with four components. These are the
real and imaginary parts of the @f$E_x@f$ and @f$E_y@f$ solution waves. With the
default setup, the output should read
 
@code
Number of active cells: 65536
Number of degrees of freedom: 263168
Output written to solution.vtk
@endcode
 
<a name="step_81-AbsorbingboundaryconditionsandthePML"></a><h3> Absorbing boundary conditions and the PML </h3>
 
 
The following images are the outputs for the imaginary @f$E_x@f$ without the
interface and with the dipole centered at @f$(0,0)@f$. In order to remove the
interface, the surface conductivity is set to 0. First, we turn off the
absorbing boundary conditions and the PML. Second, we want to see the
effect of the PML when absorbing boundary conditions apply. So we set
absorbing boundary conditions to true and leave the PML strength to 0.
Lastly, we increase the strength of the PML to 4. Change the following in
the .prm file :
 
@code
# use absorbing boundary conditions?
  set absorbing boundary condition = false
 
# position of the dipole
  set dipole position    = 0, 0
 
# strength of the PML
  set strength     = 0
 
# surface conductivity between material 1 and material 2
  set sigma              = 0, 0; 0, 0| 0, 0; 0, 0
@endcode
 
Following are the output images:
 
<table width="80%" align="center">
  <tr>
      <td align="center">
      <img src="https://dealii.org/images/steps/developer/step-81-nointerface_noabs_PML0.png" alt="Visualization of the solution of step-81 with no interface, Dirichlet boundary conditions and PML strength 0" height="210"/>
      <p> Solution with no interface, Dirichlet boundary conditions and PML strength 0.</p>
      </td>
      <td></td>
        <td align="center">
      <img src="https://dealii.org/images/steps/developer/step-81-nointerface_abs_PML0.png" alt="Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with no interface, absorbing boundary conditions and PML strength 0.</p>
            </td>
        <td></td>
            <td align="center">
      <img src="https://dealii.org/images/steps/developer/step-81-nointerface_abs_PML4.png" alt="Visualization of the solution of step-81 with no interface, absorbing boundary conditions and PML strength 4" height="210">
      <p> Solution with no interface, absorbing boundary conditions and PML strength 4.</p>
    </td>
  </tr>
</table>
 
We observe that with absorbing boundary conditions and in absence of the
PML, there is a lot of distortion and resonance (the real parts will not be
generated without a PML). This is, as we stipulated, due to reflection from
infinity. As we see, a much more coherent image is generated with an
appropriate PML.
 
<a name="step_81-SurfacePlasmonPolariton"></a><h3> Surface Plasmon Polariton </h3>
 
Now, let's generate a standing wave by adding an interface at the center.
In order to observe this effect, we offset the center of the dipole to @f$(0,
0.8)@f$ and set the surface conductivity back to @f$(0.001, 0.2)@f$:
 
@code
# position of the dipole
  set dipole position    = 0, 0.8
 
# surface conductivity between material 1 and material 2
  set sigma              = 0.001, 0.2; 0, 0| 0, 0; 0.001, 0.2
@endcode
 
Once again, we will visualize the output with absorbing boundary conditions
and PML strength 0 and with absorbing boundary conditions and PML strength
4. The following tables are the imaginary part of @f$E_x@f$ and the real part
of @f$E_x@f$.
 
<table width="80%" align="center">
  <tr>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-imagEx_noabs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
            </td>
            <td></td>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-imagEx_abs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
        </td>
        <td></td>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-imagEx_abs_PML4.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
        <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
        </td>
  </tr>
</table>
 
 
<table width="80%" align="center">
  <tr>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-realEx_noabs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
            </td>
            <td></td>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-realEx_abs_PML0.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
        </td>
        <td></td>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-realEx_abs_PML4.png" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
        <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
    </td>
  </tr>
</table>
 
The SPP is confined near the interface that we created, however without
absorbing boundary conditions, we don't observe a dissipation effect. On
adding the absorbing boundary conditions, we observe distortion and
resonance and we still don't notice any dissipation. As expected, the PML
removes the distortion and resonance. The standing wave is also dissipating
and getting absorbed within the PML, and as we increase the PML strength,
the standing wave will dissipate more within the PML ring.
 
Here are some animations to demonstrate the effect of the PML
<table width="80%" align="center">
  <tr>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-dirichlet_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
            </td>
            <td></td>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-absorbing_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
        </td>
        <td></td>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-perfectly_matched_layer_Ex.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
        <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
        </td>
  </tr>
</table>
 
 
<table width="80%" align="center">
  <tr>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-dirichlet_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
            </td>
            <td></td>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-absorbing_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 0" height="210">
        <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
        </td>
        <td></td>
        <td align="center">
        <img src="https://dealii.org/images/steps/developer/step-81-perfectly_matched_layer_Ey.gif" alt="Visualization of the solution of step-81 with an interface, absorbing boundary conditions and PML strength 4" height="210">
        <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
        </td>
  </tr>
</table>
 
<a name="step_81-Notes"></a><h3> Notes </h3>
 
 
<a name="step_81-RealandComplexMatrices"></a><h4> Real and Complex Matrices </h4>
 
As is evident from the results, we are splitting our solution matrices into
the real and the imaginary components. We started off using the @f$H^{curl}@f$
conforming Nédélec Elements, and we made two copies of the Finite Elements
in order to represent the real and the imaginary components of our input
(FE_NedelecSZ was used instead of FE_Nedelec to avoid the sign conflicts
issues present in traditional Nédélec elements). In the assembly, we create
two vectors of dimension @f$dim@f$ that assist us in extracting the real and
the imaginary components of our finite elements.
 
 
<a name="step_81-RotationsandScaling"></a><h4> Rotations and Scaling </h4>
 
As we see in our assembly, our finite element is rotated and scaled as
follows:
 
@code
const auto phi_i = real_part.value(i, q_point) - 1.0i * imag_part.value(i, q_point);
@endcode
 
This @f$\phi_i@f$ variable doesn't need to be scaled in this way, we may choose
any arbitrary scaling constants @f$a@f$ and @f$b@f$. If we choose this scaling, the
@f$\phi_j@f$ must also be modified with the same scaling, as follows:
 
@code
const auto phi_i = a*real_part.value(i, q_point) -
                                   bi * imag_part.value(i, q_point);
 
const auto phi_j = a*real_part.value(i, q_point) +
                                   bi * imag_part.value(i, q_point);
@endcode
 
Moreover, the cell_rhs need not be the real part of the rhs_value. Say if
we modify to take the imaginary part of the computed rhs_value, we must
also modify the cell_matrix accordingly to take the imaginary part of temp.
However, making these changes to both sides of the equation will not affect
our solution, and we will still be able to generate the surface plasmon
polariton.
 
@code
cell_rhs(i) += rhs_value.imag();
 
cell_matrix(i) += temp.imag();
@endcode
 
<a name="step_81-Postprocessing"></a><h4> Postprocessing </h4>
 
We will create a video demonstrating the wave in motion, which is
essentially an implementation of @f$e^{-i\omega t}(Re(E) + i*Im(E))@f$ as we
increment time. This is done by slightly changing the output function to
generate a series of .vtk files, which will represent out solution wave as
we increment time. Introduce an input variable @f$t@f$ in the output_results()
class as output_results(unsigned int t). Then change the class itself to
the following:
 
@code
template <int dim>
void Maxwell<dim>::output_results(unsigned int t)
{
  std::cout << "Running step:" << t << std::endl;
  DataOut<2> data_out;
  data_out.attach_dof_handler(dof_handler);
  Vector<double> postprocessed;
  postprocessed.reinit(solution);
  for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
    {
      if (i % 4 == 0)
        {
          postprocessed[i] = std::cos(2. * numbers::PI * 0.04 * t) * solution[i] -
                             std::sin(2. * numbers::PI * 0.04 * t) * solution[i + 1];
        }
      else if (i % 4 == 2)
        {
          postprocessed[i] = std::cos(2. * numbers::PI * 0.04 * t) * solution[i] -
                             std::sin(2. * numbers::PI * 0.04 * t) * solution[i + 1];
        }
    }
  data_out.add_data_vector(postprocessed, {"E_x", "E_y", "null0", "null1"});
  data_out.build_patches();
  const std::string filename =
    "solution-" + Utilities::int_to_string(t) + ".vtk";
  std::ofstream output(filename);
  data_out.write_vtk(output);
  std::cout << "Done running step:" << t << std::endl;
}
@endcode
 
Finally, in the run() function, replace output_results() with
@code
for (int t = 0; t <= 100; t++)
  {
    output_results(t);
  }
@endcode
 
This would generate 100 solution .vtk files, which can be opened in a group
on Paraview and then can be saved as an animation. We used FFMPEG to
generate gifs.
 
<a name="step_81-PossibilitiesforExtension"></a><h3> Possibilities for Extension </h3>
 
 
The example step could be extended in a number of different directions.
<ul>
  <li>
    The current program uses a direct solver to solve the linear system.
    This is efficient for two spatial dimensions where scattering problems
    up to a few millions degrees of freedom can be solved. In 3D, however,
    the increased stencil size of the Nedelec element pose a severe
    limiting factor on the problem size that can be computed. As an
    alternative, the idea to use iterative solvers can be entertained.
    This, however requires specialized preconditioners. For example, just
    using an iterative Krylov space solver (such as SolverGMRES) on above
    problem will requires many thousands of iterations to converge.
    Unfortunately, time-harmonic Maxwell's equations lack the usual notion
    of local smoothing properties, which renders the usual suspects, such
    as a geometric multigrid (see the Multigrid class), largely useless. A
    possible extension would be to implement an additive Schwarz preconditioner
    (based on domain decomposition, see for example
    @cite Gopalakrishnan2003), or a sweeping preconditioner (see for
    example @cite Ying2012).
  </li>
  <li>
    Another possible extension of the current program is to introduce local
    mesh refinement (either based on a residual estimator, or based on the
    dual weighted residual method, see @ref step_14 "step-14"). This is in particular of
    interest to counter the increased computational cost caused by the
    scale separation between the SPP and the dipole.
  </li>
</ul>
 *
 *
<a name="step_81-PlainProg"></a>
<h1> The plain program</h1>
@include "step-81.cc"
*/