779 * <a name=
"step_81-PerfectlyMatchedLayerClass"></a>
807 *
static_assert(dim == 2,
808 *
"The perfectly matched layer is only implemented in 2d.");
824 *
std::complex<double>
d_2,
834 *
double inner_radius;
844 *
inner_radius = 12.;
845 *
add_parameter(
"inner radius",
847 *
"inner radius of the PML shell");
849 *
add_parameter(
"outer radius",
851 *
"outer radius of the PML shell");
853 *
add_parameter(
"strength",
strength,
"strength of the PML");
858 *
typename std::complex<double>
862 *
if (radius > inner_radius)
865 *
strength * ((radius - inner_radius) * (radius - inner_radius)) /
877 *
typename std::complex<double>
881 *
if (radius > inner_radius)
883 *
const double s_bar =
885 *
((radius - inner_radius) * (radius - inner_radius) *
886 *
(radius - inner_radius)) /
889 *
return {1.0,
s_bar};
901 *
std::complex<double>
d_2,
917 *
const auto d = this->d(point);
918 *
const auto d_bar = this->d_bar(point);
928 *
const auto d = this->d(point);
929 *
const auto d_bar = this->d_bar(point);
938 *
const auto d = this->d(point);
939 *
const auto d_bar = this->d_bar(point);
948 * <a name=
"step_81-MaxwellClass"></a>
972 *
unsigned int refinements;
990 *
std::unique_ptr<FiniteElement<dim>> fe;
1002 * <a name=
"step_81-ClassTemplateDefinitionsandImplementation"></a>
1008 * <a name=
"step_81-TheConstructor"></a>
1019 * these are also known as perfectly conducting boundary conditions.
1025 * template <int dim>
1026 * Maxwell<dim>::Maxwell()
1027 * : ParameterAcceptor("Maxwell")
1028 * , dof_handler(triangulation)
1030 * ParameterAcceptor::parse_parameters_call_back.connect(
1031 * [&]() { parse_parameters_callback(); });
1034 * add_parameter("scaling", scaling, "scale of the hypercube geometry");
1037 * add_parameter("refinements",
1039 * "number of refinements of the geometry");
1042 * add_parameter("fe order", fe_order, "order of the finite element space");
1044 * quadrature_order = 1;
1045 * add_parameter("quadrature order",
1047 * "order of the quadrature");
1049 * absorbing_boundary = true;
1050 * add_parameter("absorbing boundary condition",
1051 * absorbing_boundary,
1052 * "use absorbing boundary conditions?");
1056 * template <int dim>
1057 * void Maxwell<dim>::parse_parameters_callback()
1059 * fe = std::make_unique<FESystem<dim>>(FE_NedelecSZ<dim>(fe_order), 2);
1064 * The Maxwell::make_grid() routine creates the mesh for the
1065 * computational domain which in our case is a scaled square domain.
1066 * Additionally, a material interface is introduced by setting the
1067 * material id of the upper half (@f$y>0@f$) to 1 and of the lower half
1068 * (@f$y<0@f$) of the computational domain to 2.
1069 * We are using a block decomposition into real and imaginary matrices
1070 * for the solution matrices. More details on this are available
1071 * under the Results section.
1077 * template <int dim>
1078 * void Maxwell<dim>::make_grid()
1080 * GridGenerator::hyper_cube(triangulation, -scaling, scaling);
1081 * triangulation.refine_global(refinements);
1083 * if (!absorbing_boundary)
1085 * for (auto &face : triangulation.active_face_iterators())
1086 * if (face->at_boundary())
1087 * face->set_boundary_id(1);
1090 * for (auto &cell : triangulation.active_cell_iterators())
1091 * if (cell->center()[1] > 0.)
1092 * cell->set_material_id(1);
1094 * cell->set_material_id(2);
1097 * std::cout << "Number of active cells: " << triangulation.n_active_cells()
1103 * The Maxwell::setup_system() routine follows the usual routine of
1104 * enumerating all the degrees of freedom and setting up the matrix and
1105 * vector objects to hold the system data. Enumerating is done by using
1106 * DoFHandler::distribute_dofs().
1112 * template <int dim>
1113 * void Maxwell<dim>::setup_system()
1115 * dof_handler.distribute_dofs(*fe);
1116 * std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
1119 * solution.reinit(dof_handler.n_dofs());
1120 * system_rhs.reinit(dof_handler.n_dofs());
1122 * constraints.clear();
1124 * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1126 * VectorTools::project_boundary_values_curl_conforming_l2(
1128 * 0, /* real part */
1129 * Functions::ZeroFunction<dim>(2 * dim),
1130 * 0, /* boundary id */
1132 * VectorTools::project_boundary_values_curl_conforming_l2(
1134 * dim, /* imaginary part */
1135 * Functions::ZeroFunction<dim>(2 * dim),
1136 * 0, /* boundary id */
1139 * constraints.close();
1141 * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
1142 * DoFTools::make_sparsity_pattern(dof_handler,
1145 * /* keep_constrained_dofs = */ true);
1146 * sparsity_pattern.copy_from(dsp);
1147 * system_matrix.reinit(sparsity_pattern);
1152 * This is a helper function that takes the tangential component of a tensor.
1155 * template <int dim>
1156 * DEAL_II_ALWAYS_INLINE inline Tensor<1, dim, std::complex<double>>
1157 * tangential_part(const Tensor<1, dim, std::complex<double>> &tensor,
1158 * const Tensor<1, dim> &normal)
1160 * auto result = tensor;
1161 * result[0] = normal[1] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
1162 * result[1] = -normal[0] * (tensor[0] * normal[1] - tensor[1] * normal[0]);
1169 * Assemble the stiffness matrix and the right-hand side:
1171 * A_{ij} = \int_\Omega (\mu_r^{-1}\nabla \times \varphi_j) \cdot
1172 * (\nabla\times\bar{\varphi}_i)\text{d}x
1173 * - \int_\Omega \varepsilon_r\varphi_j \cdot \bar{\varphi}_i\text{d}x
1174 * - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_j)_T) \cdot
1175 * (\bar{\varphi}_i)_T\text{do}x
1176 * - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_j)_T) \cdot
1177 * (\nabla\times(\bar{\varphi}_i)_T)\text{d}x, \f} \f{align}{
1178 * F_i = i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x - \int_\Omega
1179 * \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i}) \text{d}x.
1181 * In addition, we will be modifying the coefficients if the position of the
1182 * cell is within the PML region.
1188 * template <int dim>
1189 * void Maxwell<dim>::assemble_system()
1191 * const QGauss<dim> quadrature_formula(quadrature_order);
1192 * const QGauss<dim - 1> face_quadrature_formula(quadrature_order);
1194 * FEValues<dim, dim> fe_values(*fe,
1195 * quadrature_formula,
1196 * update_values | update_gradients |
1197 * update_quadrature_points |
1198 * update_JxW_values);
1199 * FEFaceValues<dim, dim> fe_face_values(*fe,
1200 * face_quadrature_formula,
1201 * update_values | update_gradients |
1202 * update_quadrature_points |
1203 * update_normal_vectors |
1204 * update_JxW_values);
1206 * const unsigned int dofs_per_cell = fe->dofs_per_cell;
1208 * const unsigned int n_q_points = quadrature_formula.size();
1209 * const unsigned int n_face_q_points = face_quadrature_formula.size();
1211 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1212 * Vector<double> cell_rhs(dofs_per_cell);
1213 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1217 * Next, let us assemble on the interior of the domain on the left hand
1218 * side. So we are computing
1220 * \int_\Omega (\mu_r^{-1}\nabla \times \varphi_i) \cdot
1221 * (\nabla\times\bar{\varphi}_j)\text{d}x
1223 * \int_\Omega \varepsilon_r\varphi_i \cdot \bar{\varphi}_j\text{d}x
1227 * i\int_\Omega J_a \cdot \bar{\varphi_i}\text{d}x
1228 * - \int_\Omega \mu_r^{-1} \cdot (\nabla \times \bar{\varphi_i})
1231 * In doing so, we need test functions @f$\varphi_i@f$ and @f$\varphi_j@f$, and the
1232 * curl of these test variables. We must be careful with the signs of the
1233 * imaginary parts of these complex test variables. Moreover, we have a
1234 * conditional that changes the parameters if the cell is in the PML region.
1237 * const FEValuesExtractors::Vector real_part(0);
1238 * const FEValuesExtractors::Vector imag_part(dim);
1239 * for (const auto &cell : dof_handler.active_cell_iterators())
1241 * fe_values.reinit(cell);
1246 * cell->get_dof_indices(local_dof_indices);
1247 * const auto id = cell->material_id();
1249 * const auto &quadrature_points = fe_values.get_quadrature_points();
1251 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1253 * const Point<dim> &position = quadrature_points[q_point];
1255 * auto mu_inv = parameters.mu_inv(position, id);
1256 * auto epsilon = parameters.epsilon(position, id);
1257 * const auto J_a = parameters.J_a(position, id);
1259 * const auto A = perfectly_matched_layer.a_matrix(position);
1260 * const auto B = perfectly_matched_layer.b_matrix(position);
1261 * const auto d = perfectly_matched_layer.d(position);
1263 * mu_inv = mu_inv / d;
1264 * epsilon = invert(A) * epsilon * invert(B);
1266 * for (const auto i : fe_values.dof_indices())
1268 * constexpr std::complex<double> imag{0., 1.};
1270 * const auto phi_i =
1271 * fe_values[real_part].value(i, q_point) -
1272 * imag * fe_values[imag_part].value(i, q_point);
1273 * const auto curl_phi_i =
1274 * fe_values[real_part].curl(i, q_point) -
1275 * imag * fe_values[imag_part].curl(i, q_point);
1277 * const auto rhs_value =
1278 * (imag * scalar_product(J_a, phi_i)) * fe_values.JxW(q_point);
1279 * cell_rhs(i) += rhs_value.real();
1281 * for (const auto j : fe_values.dof_indices())
1283 * const auto phi_j =
1284 * fe_values[real_part].value(j, q_point) +
1285 * imag * fe_values[imag_part].value(j, q_point);
1286 * const auto curl_phi_j =
1287 * fe_values[real_part].curl(j, q_point) +
1288 * imag * fe_values[imag_part].curl(j, q_point);
1291 * (scalar_product(mu_inv * curl_phi_j, curl_phi_i) -
1292 * scalar_product(epsilon * phi_j, phi_i)) *
1293 * fe_values.JxW(q_point);
1294 * cell_matrix(i, j) += temp.real();
1301 * Now we assemble the face and the boundary. The following loops will
1304 * - i\int_\Sigma (\sigma_r^{\Sigma}(\varphi_i)_T) \cdot
1305 * (\bar{\varphi}_j)_T\text{do}x \f} and \f{align}{
1306 * - i\int_{\partial\Omega} (\sqrt{\mu_r^{-1}\varepsilon}(\varphi_i)_T)
1307 * \cdot (\nabla\times(\bar{\varphi}_j)_T)\text{d}x,
1309 * respectively. The test variables and the PML are implemented
1310 * similarly as the domain.
1313 * const FEValuesExtractors::Vector real_part(0);
1314 * const FEValuesExtractors::Vector imag_part(dim);
1315 * for (const auto &face : cell->face_iterators())
1317 * if (face->at_boundary())
1319 * const auto id = face->boundary_id();
1322 * fe_face_values.reinit(cell, face);
1324 * for (unsigned int q_point = 0; q_point < n_face_q_points;
1327 * const auto &position = quadrature_points[q_point];
1329 * auto mu_inv = parameters.mu_inv(position, id);
1330 * auto epsilon = parameters.epsilon(position, id);
1333 * perfectly_matched_layer.a_matrix(position);
1335 * perfectly_matched_layer.b_matrix(position);
1336 * const auto d = perfectly_matched_layer.d(position);
1338 * mu_inv = mu_inv / d;
1339 * epsilon = invert(A) * epsilon * invert(B);
1341 * const auto normal =
1342 * fe_face_values.normal_vector(q_point);
1344 * for (const auto i : fe_face_values.dof_indices())
1346 * constexpr std::complex<double> imag{0., 1.};
1348 * const auto phi_i =
1349 * fe_face_values[real_part].value(i, q_point) -
1351 * fe_face_values[imag_part].value(i, q_point);
1352 * const auto phi_i_T = tangential_part(phi_i, normal);
1354 * for (const auto j : fe_face_values.dof_indices())
1356 * const auto phi_j =
1357 * fe_face_values[real_part].value(j, q_point) +
1359 * fe_face_values[imag_part].value(j, q_point);
1360 * const auto phi_j_T =
1361 * tangential_part(phi_j, normal) *
1362 * fe_face_values.JxW(q_point);
1364 * const auto prod = mu_inv * epsilon;
1365 * const auto sqrt_prod = prod;
1368 * -imag * scalar_product((sqrt_prod * phi_j_T),
1370 * cell_matrix(i, j) += temp.real();
1380 * We are on an interior face:
1383 * const auto face_index = cell->face_iterator_to_index(face);
1385 * const auto id1 = cell->material_id();
1386 * const auto id2 = cell->neighbor(face_index)->material_id();
1389 * continue; /* skip this face */
1391 * fe_face_values.reinit(cell, face);
1393 * for (unsigned int q_point = 0; q_point < n_face_q_points;
1396 * const auto &position = quadrature_points[q_point];
1398 * auto sigma = parameters.sigma(position, id1, id2);
1400 * const auto B = perfectly_matched_layer.b_matrix(position);
1401 * const auto C = perfectly_matched_layer.c_matrix(position);
1402 * sigma = invert(C) * sigma * invert(B);
1404 * const auto normal = fe_face_values.normal_vector(q_point);
1406 * for (const auto i : fe_face_values.dof_indices())
1408 * constexpr std::complex<double> imag{0., 1.};
1410 * const auto phi_i =
1411 * fe_face_values[real_part].value(i, q_point) -
1412 * imag * fe_face_values[imag_part].value(i, q_point);
1413 * const auto phi_i_T = tangential_part(phi_i, normal);
1415 * for (const auto j : fe_face_values.dof_indices())
1417 * const auto phi_j =
1418 * fe_face_values[real_part].value(j, q_point) +
1420 * fe_face_values[imag_part].value(j, q_point);
1421 * const auto phi_j_T = tangential_part(phi_j, normal);
1425 * scalar_product((sigma * phi_j_T), phi_i_T) *
1426 * fe_face_values.JxW(q_point);
1427 * cell_matrix(i, j) += temp.real();
1434 * constraints.distribute_local_to_global(
1435 * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1441 * We use a direct solver from the SparseDirectUMFPACK to solve the system
1444 * template <int dim>
1445 * void Maxwell<dim>::solve()
1447 * SparseDirectUMFPACK A_direct;
1448 * A_direct.initialize(system_matrix);
1449 * A_direct.vmult(solution, system_rhs);
1454 * The output is written into a vtk file with 4 components
1457 * template <int dim>
1458 * void Maxwell<dim>::output_results()
1460 * DataOut<2> data_out;
1461 * data_out.attach_dof_handler(dof_handler);
1462 * data_out.add_data_vector(solution,
1463 * {"real_Ex", "real_Ey", "imag_Ex", "imag_Ey"});
1464 * data_out.build_patches();
1465 * const std::string filename = "solution.vtk";
1466 * std::ofstream output(filename);
1467 * data_out.write_vtk(output);
1468 * std::cout << "Output written to " << filename << std::endl;
1472 * template <int dim>
1473 * void Maxwell<dim>::run()
1477 * assemble_system();
1482 * } // namespace Step81
1486 * The following main function calls the class @ref step_81 "step-81"(), initializes the
1487 * ParameterAcceptor, and calls the run() function.
1497 * using namespace dealii;
1499 * Step81::Maxwell<2> maxwell_2d;
1500 * ParameterAcceptor::initialize("parameters.prm");
1503 * catch (std::exception &exc)
1505 * std::cerr << std::endl
1507 * << "----------------------------------------------------"
1509 * std::cerr << "Exception on processing: " << std::endl
1510 * << exc.what() << std::endl
1511 * << "Aborting!" << std::endl
1512 * << "----------------------------------------------------"
1518 * std::cerr << std::endl
1520 * << "----------------------------------------------------"
1522 * std::cerr << "Unknown exception!" << std::endl
1523 * << "Aborting!" << std::endl
1524 * << "----------------------------------------------------"
1531<a name="step_81-Results"></a><h1>Results</h1>
1534The solution is written to a .vtk file with four components. These are the
1535real and imaginary parts of the @f$E_x@f$ and @f$E_y@f$ solution waves. With the
1536current setup, the output should read
1539Number of active cells: 4096
1540Number of degrees of freedom: 16640
1541Program ended with exit code: 0
1544<a name="step_81-AbsorbingboundaryconditionsandthePML"></a><h3> Absorbing boundary conditions and the PML </h3>
1547The following images are the outputs for the imaginary @f$E_x@f$ without the
1548interface and with the dipole centered at @f$(0,0)@f$. In order to remove the
1549interface, the surface conductivity is set to 0. First, we turn off the
1550absorbing boundary conditions and the PML. Second, we want to see the
1551effect of the PML when absorbing boundary conditions apply. So we set
1552absorbing boundary conditions to true and leave the PML strength to 0.
1553Lastly, we increase the strength of the PML to 4. Change the following in
1557# use absorbing boundary conditions?
1558 set absorbing boundary condition = false
1560# position of the dipole
1561 set dipole position = 0, 0
1563# strength of the PML
1566# surface conductivity between material 1 and material 2
1567 set sigma = 0, 0; 0, 0| 0, 0; 0, 0
1570Following are the output images:
1572<table width="80%" align="center">
1575 <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"/>
1576 <p> Solution with no interface, Dirichlet boundary conditions and PML strength 0.</p>
1580 <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">
1581 <p> Solution with no interface, absorbing boundary conditions and PML strength 0.</p>
1585 <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">
1586 <p> Solution with no interface, absorbing boundary conditions and PML strength 4.</p>
1591We observe that with absorbing boundary conditions and in absence of the
1592PML, there is a lot of distortion and resonance (the real parts will not be
1593generated without a PML). This is, as we stipulated, due to reflection from
1594infinity. As we see, a much more coherent image is generated with an
1597<a name="step_81-SurfacePlasmonPolariton"></a><h3> Surface Plasmon Polariton </h3>
1604# position of the dipole
1605 set
dipole position = 0, 0.8
1607# surface conductivity between material 1 and material 2
1608 set
sigma = 0.001, 0.2; 0, 0| 0, 0; 0.001, 0.2
1616<table width=
"80%" align=
"center">
1619 <
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">
1624 <
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">
1629 <
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">
1636<table width=
"80%" align=
"center">
1639 <
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">
1644 <
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">
1649 <
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">
1664<table width=
"80%" align=
"center">
1667 <
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">
1672 <
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">
1677 <
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">
1684<table width=
"80%" align=
"center">
1687 <
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">
1692 <
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">
1697 <
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">
1703<a name=
"step_81-Notes"></a><
h3>
Notes </
h3>
1718<a name=
"step_81-RotationsandScaling"></a><
h4> Rotations
and Scaling </
h4>
1728any arbitrary scaling constants @f$a@f$ and @f$b@f$. If we choose this scaling, the
1729@f$\phi_j@f$ must also be modified with the same scaling, as follows:
1732const auto phi_i = a*real_part.value(i, q_point) -
1733 bi * imag_part.value(i, q_point);
1735const auto phi_j = a*real_part.value(i, q_point) +
1736 bi * imag_part.value(i, q_point);
1739Moreover, the cell_rhs need not be the real part of the rhs_value. Say if
1740we modify to take the imaginary part of the computed rhs_value, we must
1741also modify the cell_matrix accordingly to take the imaginary part of temp.
1742However, making these changes to both sides of the equation will not affect
1743our solution, and we will still be able to generate the surface plasmon
1747cell_rhs(i) += rhs_value.imag();
1749cell_matrix(i) += temp.imag();
1752<a name="step_81-Postprocessing"></a><h4> Postprocessing </h4>
1754We will create a video demonstrating the wave in motion, which is
1755essentially an implementation of @f$e^{-i\omega t}(Re(E) + i*Im(E))@f$ as we
1756increment time. This is done by slightly changing the output function to
1757generate a series of .vtk files, which will represent out solution wave as
1758we increment time. Introduce an input variable @f$t@f$ in the output_results()
1759class as output_results(unsigned int t). Then change the class itself to
1764void Maxwell<dim>::output_results(unsigned int t)
1766 std::cout << "Running step:" << t << std::endl;
1767 DataOut<2> data_out;
1768 data_out.attach_dof_handler(dof_handler);
1769 Vector<double> postprocessed;
1770 postprocessed.reinit(solution);
1771 for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1775 postprocessed[i] = std::cos(2 * M_PI * 0.04 * t) * solution[i] -
1776 std::sin(2 * M_PI * 0.04 * t) * solution[i + 1];
1778 else if (i % 4 == 2)
1780 postprocessed[i] = std::cos(2 * M_PI * 0.04 * t) * solution[i] -
1781 std::sin(2 * M_PI * 0.04 * t) * solution[i + 1];
1784 data_out.add_data_vector(postprocessed, {"E_x", "E_y", "null0", "null1"});
1785 data_out.build_patches();
1786 const std::string filename =
1787 "solution-" + Utilities::int_to_string(t) + ".vtk";
1788 std::ofstream output(filename);
1789 data_out.write_vtk(output);
1790 std::cout << "Done running step:" << t << std::endl;
1794Finally, in the run() function, replace output_results() with
1796for (int t = 0; t <= 100; t++)
1802This would generate 100 solution .vtk files, which can be opened in a group
1803on Paraview and then can be saved as an animation. We used FFMPEG to
1806<a name="step_81-PossibilitiesforExtension"></a><h3> Possibilities for Extension </h3>
1809The example step could be extended in a number of different directions.
1812 The current program uses a direct solver to solve the linear system.
1813 This is efficient for two spatial dimensions where scattering problems
1814 up to a few millions degrees of freedom can be solved. In 3D, however,
1815 the increased stencil size of the Nedelec element pose a severe
1816 limiting factor on the problem size that can be computed. As an
1817 alternative, the idea to use iterative solvers can be entertained.
1818 This, however requires specialized preconditioners. For example, just
1819 using an iterative Krylov space solver (such as SolverGMRES) on above
1820 problem will requires many thousands of iterations to converge.
1839<a name=
"step_81-PlainProg"></a>
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
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
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 > &)