741 *
return (material == 1 ? epsilon_1 : epsilon_2);
745 * std::complex<double> Parameters<dim>::mu_inv(
const Point<dim> & ,
748 *
return (material == 1 ? mu_inv_1 : mu_inv_2);
752 *
typename Parameters<dim>::rank2_type
757 *
return (left == right ? rank2_type() : sigma_tensor);
761 *
typename Parameters<dim>::rank1_type
765 *
const auto distance = (dipole_position -
point).
norm() / dipole_radius;
769 *
std::cos(distance * M_PI / 2.) / (M_PI / 2. - 2. / M_PI) /
770 * dipole_radius / dipole_radius;
771 * J_a = dipole_strength * dipole_orientation *
scale;
778 * <a name=
"step_81-PerfectlyMatchedLayerClass"></a>
779 * <h4>PerfectlyMatchedLayer Class</h4>
781 * implements the transformation matrices used to modify the permittivity
782 * and permeability tensors supplied from the Parameters
class. The
783 * actual transformation of the material tensors will be done in the
784 * assembly
loop. The radii and the strength of the PML is specified, and
785 * the coefficients will be modified
using transformation matrices within
786 * the PML region. The radii and strength of the PML are editable through
787 * a .prm file. The rotation function @f$T_{exer}@f$ is the same as
788 * introduced in the perfectly matched layer section of the introduction.
789 * Similarly, the matrices A, B and
C are defined as follows
791 * A = T_{e_xe_r}^{-1}
792 * \text{diag}\left(\frac{1}{\bar{
d}^2},\frac{1}{d\bar{
d}}\right)T_{e_xe_r},\qquad
793 * B = T_{e_xe_r}^{-1} \text{diag}\left(d,\bar{
d}\right)T_{e_xe_r},\qquad
794 *
C = T_{e_xe_r}^{-1} \text{diag}\left(\frac{1}{\bar{
d}},\frac{1}{
d}\right)
806 *
static_assert(dim == 2,
807 *
"The perfectly matched layer is only implemented in 2d.");
809 * Parameters<dim> parameters;
815 * PerfectlyMatchedLayer();
819 * std::complex<double> d_bar(
const Point<dim> point);
822 * rank2_type rotation(std::complex<double> d_1,
823 * std::complex<double> d_2,
826 * rank2_type a_matrix(
const Point<dim> point);
828 * rank2_type b_matrix(
const Point<dim> point);
830 * rank2_type c_matrix(
const Point<dim> point);
833 *
double inner_radius;
834 *
double outer_radius;
840 * PerfectlyMatchedLayer<dim>::PerfectlyMatchedLayer()
843 * inner_radius = 12.;
844 * add_parameter(
"inner radius",
846 *
"inner radius of the PML shell");
847 * outer_radius = 20.;
848 * add_parameter(
"outer radius",
850 *
"outer radius of the PML shell");
852 * add_parameter(
"strength", strength,
"strength of the PML");
857 *
typename std::complex<double>
858 * PerfectlyMatchedLayer<dim>::d(
const Point<dim> point)
861 *
if (radius > inner_radius)
864 * strength * ((radius - inner_radius) * (radius - inner_radius)) /
865 * ((outer_radius - inner_radius) * (outer_radius - inner_radius));
876 *
typename std::complex<double>
877 * PerfectlyMatchedLayer<dim>::d_bar(
const Point<dim> point)
880 *
if (radius > inner_radius)
882 *
const double s_bar =
884 * ((radius - inner_radius) * (radius - inner_radius) *
885 * (radius - inner_radius)) /
886 * (radius * (outer_radius - inner_radius) *
887 * (outer_radius - inner_radius));
888 *
return {1.0, s_bar};
898 *
typename PerfectlyMatchedLayer<dim>::rank2_type
899 * PerfectlyMatchedLayer<dim>::rotation(std::complex<double> d_1,
900 * std::complex<double> d_2,
905 * result[0][1] =
point[0] *
point[1] * (d_1 - d_2);
906 * result[1][0] =
point[0] *
point[1] * (d_1 - d_2);
913 *
typename PerfectlyMatchedLayer<dim>::rank2_type
914 * PerfectlyMatchedLayer<dim>::a_matrix(
const Point<dim> point)
916 *
const auto d = this->d(point);
917 *
const auto d_bar = this->d_bar(point);
918 *
return invert(rotation(d * d, d * d_bar, point)) *
919 * rotation(d * d, d * d_bar, point);
924 *
typename PerfectlyMatchedLayer<dim>::rank2_type
925 * PerfectlyMatchedLayer<dim>::b_matrix(
const Point<dim> point)
927 *
const auto d = this->d(point);
928 *
const auto d_bar = this->d_bar(point);
929 *
return invert(rotation(d, d_bar, point)) * rotation(d, d_bar, point);
934 *
typename PerfectlyMatchedLayer<dim>::rank2_type
935 * PerfectlyMatchedLayer<dim>::c_matrix(
const Point<dim> point)
937 *
const auto d = this->d(point);
938 *
const auto d_bar = this->d_bar(point);
939 *
return invert(rotation(1. / d_bar, 1. / d, point)) *
940 * rotation(1. / d_bar, 1. / d, point);
947 * <a name=
"step_81-MaxwellClass"></a>
948 * <h4>Maxwell Class</h4>
949 * At
this point we are ready to declare all the major building blocks of
950 * the finite element program which consists of the usual setup and
951 * assembly routines. Most of the structure has already been introduced
952 * in previous tutorial programs. The Maxwell
class also holds private
953 * instances of the Parameters and PerfectlyMatchedLayers classes
954 * introduced above. The
default values of these parameters are set to
955 * show us a standing wave with absorbing boundary conditions and a PML.
971 *
unsigned int refinements;
972 *
unsigned int fe_order;
973 *
unsigned int quadrature_order;
974 *
bool absorbing_boundary;
976 *
void parse_parameters_callback();
978 *
void setup_system();
979 *
void assemble_system();
981 *
void output_results();
983 * Parameters<dim> parameters;
984 * PerfectlyMatchedLayer<dim> perfectly_matched_layer;
989 * std::unique_ptr<FiniteElement<dim>> fe;
1001 * <a name=
"step_81-ClassTemplateDefinitionsandImplementation"></a>
1002 * <h3>Class Template Definitions and Implementation</h3>
1007 * <a name=
"step_81-TheConstructor"></a>
1008 * <h4>The Constructor</h4>
1009 * The Constructor simply consists of
default initialization a number of
1010 * discretization parameters (such as the domain size, mesh refinement,
1011 * and the order of finite elements and quadrature) and declaring a
1013 * these can be modified by editing the .prm file. Absorbing boundary
1014 * conditions can be controlled with the absorbing_boundary
boolean. If
1015 * absorbing boundary conditions are disabled we simply enforce
1016 * homogeneous Dirichlet conditions on the tangential component of the
1017 * electric field. In the context of time-harmonic Maxwell
's equations
1018 * these are also known as perfectly conducting boundary conditions.
1024 * template <int dim>
1025 * Maxwell<dim>::Maxwell()
1026 * : ParameterAcceptor("Maxwell")
1027 * , dof_handler(triangulation)
1029 * ParameterAcceptor::parse_parameters_call_back.connect(
1030 * [&]() { parse_parameters_callback(); });
1033 * add_parameter("scaling", scaling, "scale of the hypercube geometry");
1036 * add_parameter("refinements",
1038 * "number of refinements of the geometry");
1041 * add_parameter("fe order", fe_order, "order of the finite element space");
1043 * quadrature_order = 1;
1044 * add_parameter("quadrature order",
1046 * "order of the quadrature");
1048 * absorbing_boundary = true;
1049 * add_parameter("absorbing boundary condition",
1050 * absorbing_boundary,
1051 * "use absorbing boundary conditions?");
1055 * template <int dim>
1056 * void Maxwell<dim>::parse_parameters_callback()
1058 * fe = std::make_unique<FESystem<dim>>(FE_NedelecSZ<dim>(fe_order), 2);
1063 * The Maxwell::make_grid() routine creates the mesh for the
1064 * computational domain which in our case is a scaled square domain.
1065 * Additionally, a material interface is introduced by setting the
1066 * material id of the upper half (@f$y>0@f$) to 1 and of the lower half
1067 * (@f$y<0@f$) of the computational domain to 2.
1068 * We are using a block decomposition into real and imaginary matrices
1069 * for the solution matrices. More details on this are available
1070 * under the Results section.
1076 * template <int dim>
1077 * void Maxwell<dim>::make_grid()
1079 * GridGenerator::hyper_cube(triangulation, -scaling, scaling);
1080 * triangulation.refine_global(refinements);
1082 * if (!absorbing_boundary)
1084 * for (auto &face : triangulation.active_face_iterators())
1085 * if (face->at_boundary())
1086 * face->set_boundary_id(1);
1089 * for (auto &cell : triangulation.active_cell_iterators())
1090 * if (cell->center()[1] > 0.)
1091 * cell->set_material_id(1);
1093 * cell->set_material_id(2);
1096 * std::cout << "Number of active cells: " << triangulation.n_active_cells()
1102 * The Maxwell::setup_system() routine follows the usual routine of
1103 * enumerating all the degrees of freedom and setting up the matrix and
1104 * vector objects to hold the system data. Enumerating is done by using
1105 * DoFHandler::distribute_dofs().
1111 * template <int dim>
1112 * void Maxwell<dim>::setup_system()
1114 * dof_handler.distribute_dofs(*fe);
1115 * std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
1118 * solution.reinit(dof_handler.n_dofs());
1119 * system_rhs.reinit(dof_handler.n_dofs());
1121 * constraints.clear();
1123 * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1125 * VectorTools::project_boundary_values_curl_conforming_l2(
1127 * 0, /* real part */
1128 * Functions::ZeroFunction<dim>(2 * dim),
1129 * 0, /* boundary id */
1131 * VectorTools::project_boundary_values_curl_conforming_l2(
1133 * dim, /* imaginary part */
1134 * Functions::ZeroFunction<dim>(2 * dim),
1135 * 0, /* boundary id */
1138 * constraints.close();
1140 * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
1141 * DoFTools::make_sparsity_pattern(dof_handler,
1144 * /* keep_constrained_dofs = */ true);
1145 * sparsity_pattern.copy_from(dsp);
1146 * system_matrix.reinit(sparsity_pattern);
1151 * This is a helper function that takes the tangential component of a tensor.
1154 * template <int dim>
1155 * DEAL_II_ALWAYS_INLINE inline Tensor<1, dim, std::complex<double>>
1156 * tangential_part(const Tensor<1, dim, std::complex<double>> &tensor,
1157 * const Tensor<1, dim> &normal)
1159 * auto result = tensor;
1160 * result[0] = normal[1] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
1161 * result[1] = -normal[0] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
1168 * Assemble the stiffness matrix and the right-hand side:
1170 * A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_j) \cdot
1171 * (\nabla\times\bar{\varphi}_i)\text{d}x
1172 * - \int_\Omega \varepsilon_r\varphi_j \cdot \bar{\varphi}_i\text{d}x
1173 * - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_j)_T) \cdot
1174 * (\bar{\varphi}_i)_T\text{do}x
1175 * - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_j)_T) \cdot
1176 * (\nabla\times(\bar{\varphi}_i)_T)\text{d}x, \f} \f{align}{
1177 * F_i = i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x - \int_\Omega
1178 * \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x.
1180 * In addition, we will be modifying the coefficients if the position of the
1181 * cell is within the PML region.
1187 * template <int dim>
1188 * void Maxwell<dim>::assemble_system()
1190 * const QGauss<dim> quadrature_formula(quadrature_order);
1191 * const QGauss<dim - 1> face_quadrature_formula(quadrature_order);
1193 * FEValues<dim, dim> fe_values(*fe,
1194 * quadrature_formula,
1195 * update_values | update_gradients |
1196 * update_quadrature_points |
1197 * update_JxW_values);
1198 * FEFaceValues<dim, dim> fe_face_values(*fe,
1199 * face_quadrature_formula,
1200 * update_values | update_gradients |
1201 * update_quadrature_points |
1202 * update_normal_vectors |
1203 * update_JxW_values);
1205 * const unsigned int dofs_per_cell = fe->dofs_per_cell;
1207 * const unsigned int n_q_points = quadrature_formula.size();
1208 * const unsigned int n_face_q_points = face_quadrature_formula.size();
1210 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1211 * Vector<double> cell_rhs(dofs_per_cell);
1212 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1216 * Next, let us assemble on the interior of the domain on the left hand
1217 * side. So we are computing
1219 * \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot
1220 * (\nabla\times\bar{\varphi}_j)\text{d}x
1222 * \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x
1226 * i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x
1227 * - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i})
1230 * In doing so, we need test functions @f$\varphi_i@f$ and @f$\varphi_j@f$, and the
1231 * curl of these test variables. We must be careful with the signs of the
1232 * imaginary parts of these complex test variables. Moreover, we have a
1233 * conditional that changes the parameters if the cell is in the PML region.
1236 * const FEValuesExtractors::Vector real_part(0);
1237 * const FEValuesExtractors::Vector imag_part(dim);
1238 * for (const auto &cell : dof_handler.active_cell_iterators())
1240 * fe_values.reinit(cell);
1245 * cell->get_dof_indices(local_dof_indices);
1246 * const auto id = cell->material_id();
1248 * const auto &quadrature_points = fe_values.get_quadrature_points();
1250 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1252 * const Point<dim> &position = quadrature_points[q_point];
1254 * auto mu_inv = parameters.mu_inv(position, id);
1255 * auto epsilon = parameters.epsilon(position, id);
1256 * const auto J_a = parameters.J_a(position, id);
1258 * const auto A = perfectly_matched_layer.a_matrix(position);
1259 * const auto B = perfectly_matched_layer.b_matrix(position);
1260 * const auto d = perfectly_matched_layer.d(position);
1262 * mu_inv = mu_inv / d;
1263 * epsilon = invert(A) * epsilon * invert(B);
1265 * for (const auto i : fe_values.dof_indices())
1267 * constexpr std::complex<double> imag{0., 1.};
1269 * const auto phi_i =
1270 * fe_values[real_part].value(i, q_point) -
1271 * imag * fe_values[imag_part].value(i, q_point);
1272 * const auto curl_phi_i =
1273 * fe_values[real_part].curl(i, q_point) -
1274 * imag * fe_values[imag_part].curl(i, q_point);
1276 * const auto rhs_value =
1277 * (imag * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point);
1278 * cell_rhs(i) += rhs_value.real();
1280 * for (const auto j : fe_values.dof_indices())
1282 * const auto phi_j =
1283 * fe_values[real_part].value(j, q_point) +
1284 * imag * fe_values[imag_part].value(j, q_point);
1285 * const auto curl_phi_j =
1286 * fe_values[real_part].curl(j, q_point) +
1287 * imag * fe_values[imag_part].curl(j, q_point);
1290 * (scalar_product(mu_inv * curl_phi_j, curl_phi_i) -
1291 * scalar_product(epsilon * phi_j, phi_i)) *
1292 * fe_values.JxW(q_point);
1293 * cell_matrix(i, j) += temp.real();
1300 * Now we assemble the face and the boundary. The following loops will
1303 * - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot
1304 * (\bar{\varphi}_j)_T\text{do}x \f} and \f{align}{
1305 * - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T)
1306 * \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x,
1308 * respectively. The test variables and the PML are implemented
1309 * similarly as the domain.
1312 * const FEValuesExtractors::Vector real_part(0);
1313 * const FEValuesExtractors::Vector imag_part(dim);
1314 * for (const auto &face : cell->face_iterators())
1316 * if (face->at_boundary())
1318 * const auto id = face->boundary_id();
1321 * fe_face_values.reinit(cell, face);
1323 * for (unsigned int q_point = 0; q_point < n_face_q_points;
1326 * const auto &position = quadrature_points[q_point];
1328 * auto mu_inv = parameters.mu_inv(position, id);
1329 * auto epsilon = parameters.epsilon(position, id);
1332 * perfectly_matched_layer.a_matrix(position);
1334 * perfectly_matched_layer.b_matrix(position);
1335 * const auto d = perfectly_matched_layer.d(position);
1337 * mu_inv = mu_inv / d;
1338 * epsilon = invert(A) * epsilon * invert(B);
1340 * const auto normal =
1341 * fe_face_values.normal_vector(q_point);
1343 * for (const auto i : fe_face_values.dof_indices())
1345 * constexpr std::complex<double> imag{0., 1.};
1347 * const auto phi_i =
1348 * fe_face_values[real_part].value(i, q_point) -
1350 * fe_face_values[imag_part].value(i, q_point);
1351 * const auto phi_i_T = tangential_part(phi_i, normal);
1353 * for (const auto j : fe_face_values.dof_indices())
1355 * const auto phi_j =
1356 * fe_face_values[real_part].value(j, q_point) +
1358 * fe_face_values[imag_part].value(j, q_point);
1359 * const auto phi_j_T =
1360 * tangential_part(phi_j, normal) *
1361 * fe_face_values.JxW(q_point);
1363 * const auto prod = mu_inv * epsilon;
1364 * const auto sqrt_prod = prod;
1367 * -imag * scalar_product((sqrt_prod * phi_j_T),
1369 * cell_matrix(i, j) += temp.real();
1379 * We are on an interior face:
1382 * const auto face_index = cell->face_iterator_to_index(face);
1384 * const auto id1 = cell->material_id();
1385 * const auto id2 = cell->neighbor(face_index)->material_id();
1388 * continue; /* skip this face */
1390 * fe_face_values.reinit(cell, face);
1392 * for (unsigned int q_point = 0; q_point < n_face_q_points;
1395 * const auto &position = quadrature_points[q_point];
1397 * auto sigma = parameters.sigma(position, id1, id2);
1399 * const auto B = perfectly_matched_layer.b_matrix(position);
1400 * const auto C = perfectly_matched_layer.c_matrix(position);
1401 * sigma = invert(C) * sigma * invert(B);
1403 * const auto normal = fe_face_values.normal_vector(q_point);
1405 * for (const auto i : fe_face_values.dof_indices())
1407 * constexpr std::complex<double> imag{0., 1.};
1409 * const auto phi_i =
1410 * fe_face_values[real_part].value(i, q_point) -
1411 * imag * fe_face_values[imag_part].value(i, q_point);
1412 * const auto phi_i_T = tangential_part(phi_i, normal);
1414 * for (const auto j : fe_face_values.dof_indices())
1416 * const auto phi_j =
1417 * fe_face_values[real_part].value(j, q_point) +
1419 * fe_face_values[imag_part].value(j, q_point);
1420 * const auto phi_j_T = tangential_part(phi_j, normal);
1424 * scalar_product((sigma * phi_j_T), phi_i_T) *
1425 * fe_face_values.JxW(q_point);
1426 * cell_matrix(i, j) += temp.real();
1433 * constraints.distribute_local_to_global(
1434 * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1440 * We use a direct solver from the SparseDirectUMFPACK to solve the system
1443 * template <int dim>
1444 * void Maxwell<dim>::solve()
1446 * SparseDirectUMFPACK A_direct;
1447 * A_direct.initialize(system_matrix);
1448 * A_direct.vmult(solution, system_rhs);
1453 * The output is written into a vtk file with 4 components
1456 * template <int dim>
1457 * void Maxwell<dim>::output_results()
1459 * DataOut<2> data_out;
1460 * data_out.attach_dof_handler(dof_handler);
1461 * data_out.add_data_vector(solution,
1462 * {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"});
1463 * data_out.build_patches();
1464 * const std::string filename = "solution.vtk";
1465 * std::ofstream output(filename);
1466 * data_out.write_vtk(output);
1467 * std::cout << "Output written to " << filename << std::endl;
1471 * template <int dim>
1472 * void Maxwell<dim>::run()
1476 * assemble_system();
1481 * } // namespace Step81
1485 * The following main function calls the class @ref step_81 "step-81"(), initializes the
1486 * ParameterAcceptor, and calls the run() function.
1496 * using namespace dealii;
1498 * Step81::Maxwell<2> maxwell_2d;
1499 * ParameterAcceptor::initialize("parameters.prm");
1502 * catch (std::exception &exc)
1504 * std::cerr << std::endl
1506 * << "----------------------------------------------------"
1508 * std::cerr << "Exception on processing: " << std::endl
1509 * << exc.what() << std::endl
1510 * << "Aborting!" << std::endl
1511 * << "----------------------------------------------------"
1517 * std::cerr << std::endl
1519 * << "----------------------------------------------------"
1521 * std::cerr << "Unknown exception!" << std::endl
1522 * << "Aborting!" << std::endl
1523 * << "----------------------------------------------------"
1530<a name="step_81-Results"></a><h1>Results</h1>
1533The solution is written to a .vtk file with four components. These are the
1534real and imaginary parts of the @f$E_x@f$ and @f$E_y@f$ solution waves. With the
1535current setup, the output should read
1538Number of active cells: 4096
1539Number of degrees of freedom: 16640
1540Program ended with exit code: 0
1543<a name="step_81-AbsorbingboundaryconditionsandthePML"></a><h3> Absorbing boundary conditions and the PML </h3>
1546The following images are the outputs for the imaginary @f$E_x@f$ without the
1547interface and with the dipole centered at @f$(0,0)@f$. In order to remove the
1548interface, the surface conductivity is set to 0. First, we turn off the
1549absorbing boundary conditions and the PML. Second, we want to see the
1550effect of the PML when absorbing boundary conditions apply. So we set
1551absorbing boundary conditions to true and leave the PML strength to 0.
1552Lastly, we increase the strength of the PML to 4. Change the following in
1556# use absorbing boundary conditions?
1557 set absorbing boundary condition = false
1559# position of the dipole
1560 set dipole position = 0, 0
1562# strength of the PML
1565# surface conductivity between material 1 and material 2
1566 set sigma = 0, 0; 0, 0| 0, 0; 0, 0
1569Following are the output images:
1571<table width="80%" align="center">
1574 <img src="https://www.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"/>
1575 <p> Solution with no interface, Dirichlet boundary conditions and PML strength 0.</p>
1579 <img src="https://www.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">
1580 <p> Solution with no interface, absorbing boundary conditions and PML strength 0.</p>
1584 <img src="https://www.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">
1585 <p> Solution with no interface, absorbing boundary conditions and PML strength 4.</p>
1590We observe that with absorbing boundary conditions and in absence of the
1591PML, there is a lot of distortion and resonance (the real parts will not be
1592generated without a PML). This is, as we stipulated, due to reflection from
1593infinity. As we see, a much more coherent image is generated with an
1596<a name="step_81-SurfacePlasmonPolariton"></a><h3> Surface Plasmon Polariton </h3>
1598Now, let's generate a standing wave by adding an
interface at the
center.
1599In order to observe this effect, we offset the
center of the dipole to @f$(0,
16000.8)@f$ and set the surface conductivity back to @f$(0.001, 0.2)@f$:
1603# position of the dipole
1604 set dipole position = 0, 0.8
1606# surface conductivity between material 1 and material 2
1607 set sigma = 0.001, 0.2; 0, 0| 0, 0; 0.001, 0.2
1610Once again, we will visualize the output with absorbing boundary conditions
1611and PML strength 0 and with absorbing boundary conditions and PML strength
16124. The following tables are the imaginary part of @f$E_x@f$ and the real part
1615<table width=
"80%" align=
"center">
1618 <img src=
"https://www.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">
1619 <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
1623 <img src=
"https://www.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">
1624 <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
1628 <img src=
"https://www.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">
1629 <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
1635<table width=
"80%" align=
"center">
1638 <img src=
"https://www.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">
1639 <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
1643 <img src=
"https://www.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">
1644 <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
1648 <img src=
"https://www.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">
1649 <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
1654The SPP is confined near the
interface that we created, however without
1655absorbing boundary conditions, we don't observe a dissipation effect. On
1656adding the absorbing boundary conditions, we observe distortion and
1657resonance and we still don't notice any dissipation. As expected, the PML
1658removes the distortion and resonance. The standing wave is also dissipating
1659and getting absorbed within the PML, and as we increase the PML strength,
1660the standing wave will dissipate more within the PML ring.
1662Here are some animations to demonstrate the effect of the PML
1663<table width=
"80%" align=
"center">
1666 <img src=
"https://www.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">
1667 <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
1671 <img src=
"https://www.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">
1672 <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
1676 <img src=
"https://www.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">
1677 <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
1683<table width=
"80%" align=
"center">
1686 <img src=
"https://www.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">
1687 <p> Solution with an interface, Dirichlet boundary conditions and PML strength 0.</p>
1691 <img src=
"https://www.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">
1692 <p> Solution with an interface, absorbing boundary conditions and PML strength 0.</p>
1696 <img src=
"https://www.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">
1697 <p> Solution with an interface, absorbing boundary conditions and PML strength 4.</p>
1702<a name=
"step_81-Notes"></a><h3> Notes </h3>
1705<a name=
"step_81-RealandComplexMatrices"></a><h4> Real and Complex Matrices </h4>
1707As is evident from the results, we are splitting our solution matrices into
1708the real and the imaginary components. We started off
using the @f$H^{curl}@f$
1709conforming Nédélec Elements, and we made two copies of the Finite Elements
1710in order to represent the real and the imaginary components of our input
1712issues present in traditional Nédélec elements). In the assembly, we create
1713two vectors of dimension @f$dim@f$ that assist us in extracting the real and
1714the imaginary components of our finite elements.
1717<a name=
"step_81-RotationsandScaling"></a><h4> Rotations and Scaling </h4>
1719As we see in our assembly, our finite element is rotated and scaled as
1723const auto phi_i =
real_part.value(i, q_point) - 1.0i * imag_part.value(i, q_point);
1726This @f$\phi_i@f$ variable doesn
't need to be scaled in this way, we may choose
1727any arbitrary scaling constants @f$a@f$ and @f$b@f$. If we choose this scaling, the
1728@f$\phi_j@f$ must also be modified with the same scaling, as follows:
1731const auto phi_i = a*real_part.value(i, q_point) -
1732 bi * imag_part.value(i, q_point);
1734const auto phi_j = a*real_part.value(i, q_point) +
1735 bi * imag_part.value(i, q_point);
1738Moreover, the cell_rhs need not be the real part of the rhs_value. Say if
1739we modify to take the imaginary part of the computed rhs_value, we must
1740also modify the cell_matrix accordingly to take the imaginary part of temp.
1741However, making these changes to both sides of the equation will not affect
1742our solution, and we will still be able to generate the surface plasmon
1746cell_rhs(i) += rhs_value.imag();
1748cell_matrix(i) += temp.imag();
1751<a name="step_81-Postprocessing"></a><h4> Postprocessing </h4>
1753We will create a video demonstrating the wave in motion, which is
1754essentially an implementation of @f$e^{-i\omega t}(Re(E) + i*Im(E))@f$ as we
1755increment time. This is done by slightly changing the output function to
1756generate a series of .vtk files, which will represent out solution wave as
1757we increment time. Introduce an input variable @f$t@f$ in the output_results()
1758class as output_results(unsigned int t). Then change the class itself to
1763void Maxwell<dim>::output_results(unsigned int t)
1765 std::cout << "Running step:" << t << std::endl;
1766 DataOut<2> data_out;
1767 data_out.attach_dof_handler(dof_handler);
1768 Vector<double> postprocessed;
1769 postprocessed.reinit(solution);
1770 for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1774 postprocessed[i] = std::cos(2 * M_PI * 0.04 * t) * solution[i] -
1775 std::sin(2 * M_PI * 0.04 * t) * solution[i + 1];
1777 else if (i % 4 == 2)
1779 postprocessed[i] = std::cos(2 * M_PI * 0.04 * t) * solution[i] -
1780 std::sin(2 * M_PI * 0.04 * t) * solution[i + 1];
1783 data_out.add_data_vector(postprocessed, {"E_x", "E_y", "null0", "null1"});
1784 data_out.build_patches();
1785 const std::string filename =
1786 "solution-" + Utilities::int_to_string(t) + ".vtk";
1787 std::ofstream output(filename);
1788 data_out.write_vtk(output);
1789 std::cout << "Done running step:" << t << std::endl;
1793Finally, in the run() function, replace output_results() with
1795for (int t = 0; t <= 100; t++)
1801This would generate 100 solution .vtk files, which can be opened in a group
1802on Paraview and then can be saved as an animation. We used FFMPEG to
1805<a name="step_81-PossibilitiesforExtension"></a><h3> Possibilities for Extension </h3>
1808The example step could be extended in a number of different directions.
1811 The current program uses a direct solver to solve the linear system.
1812 This is efficient for two spatial dimensions where scattering problems
1813 up to a few millions degrees of freedom can be solved. In 3D, however,
1814 the increased stencil size of the Nedelec element pose a severe
1815 limiting factor on the problem size that can be computed. As an
1816 alternative, the idea to use iterative solvers can be entertained.
1817 This, however requires specialized preconditioners. For example, just
1818 using an iterative Krylov space solver (such as SolverGMRES) on above
1819 problem will requires many thousands of iterations to converge.
1820 Unfortunately, time-harmonic Maxwell's equations lack the usual notion
1821 of local smoothing properties, which renders the usual suspects, such
1822 as a geometric multigrid (see the
Multigrid class), largely useless. A
1823 possible extension would be to implement an additive Schwarz preconditioner
1824 (based on domain decomposition, see
for example
1825 @cite Gopalakrishnan2003), or a sweeping preconditioner (see
for
1826 example @cite Ying2012).
1829 Another possible extension of the current program is to introduce local
1830 mesh refinement (either based on a residual estimator, or based on the
1831 dual weighted residual method, see @ref step_14
"step-14"). This is in particular of
1832 interest to counter the increased computational cost caused by the
1833 scale separation between the SPP and the dipole.
1838<a name=
"step_81-PlainProg"></a>
1839<h1> The plain program</h1>
1840@include
"step-81.cc"
void add_parameter(const std::string &entry, ParameterType ¶meter, const std::string &documentation="", ParameterHandler &prm_=prm, const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern())
numbers::NumberTraits< Number >::real_type norm() const
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)