Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-62.h
Go to the documentation of this file.
1 ) const
1016  * {
1017  * @endcode
1018  *
1019  * The speed of sound is defined by
1020  * @f[
1021  * c = \frac{K_e}{\rho}
1022  * @f]
1023  * where @f$K_e@f$ is the effective elastic constant and @f$\rho@f$ the density.
1024  * Here we consider the case in which the waveguide width is much smaller
1025  * than the wavelength. In this case it can be shown that for the two
1026  * dimensional case
1027  * @f[
1028  * K_e = 4\mu\frac{\lambda +\mu}{\lambda+2\mu}
1029  * @f]
1030  * and for the three dimensional case @f$K_e@f$ is equal to the Young's modulus.
1031  * @f[
1032  * K_e = \mu\frac{3\lambda +2\mu}{\lambda+\mu}
1033  * @f]
1034  *
1035  * @code
1036  * double elastic_constant;
1037  * if (dim == 2)
1038  * {
1039  * elastic_constant = 4 * mu * (lambda + mu) / (lambda + 2 * mu);
1040  * }
1041  * else if (dim == 3)
1042  * {
1043  * elastic_constant = mu * (3 * lambda + 2 * mu) / (lambda + mu);
1044  * }
1045  * else
1046  * {
1047  * Assert(false, ExcInternalError());
1048  * }
1049  * const double material_a_speed_of_sound =
1050  * std::sqrt(elastic_constant / material_a_rho);
1051  * const double material_a_wavelength =
1052  * material_a_speed_of_sound / cavity_resonance_frequency;
1053  * const double material_b_speed_of_sound =
1054  * std::sqrt(elastic_constant / material_b_rho);
1055  * const double material_b_wavelength =
1056  * material_b_speed_of_sound / cavity_resonance_frequency;
1057  *
1058  * @endcode
1059  *
1060  * The density @f$\rho@f$ takes the following form
1061  * <img alt="Phononic superlattice cavity"
1062  * src="https://www.dealii.org/images/steps/developer/step-62.04.svg"
1063  * height="200" />
1064  * where the brown color represents material_a and the green color
1065  * represents material_b.
1066  *
1067  * @code
1068  * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
1069  * {
1070  * const double layer_transition_center =
1071  * material_a_wavelength / 2 +
1072  * idx * (material_b_wavelength / 4 + material_a_wavelength / 4);
1073  * if (std::abs(p[0]) >=
1074  * (layer_transition_center - average_rho_width / 2) &&
1075  * std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1076  * {
1077  * const double coefficient =
1078  * (std::abs(p[0]) -
1079  * (layer_transition_center - average_rho_width / 2)) /
1080  * average_rho_width;
1081  * return (1 - coefficient) * material_a_rho +
1082  * coefficient * material_b_rho;
1083  * }
1084  * }
1085  *
1086  * @endcode
1087  *
1088  * Here we define the
1089  * [subpixel
1090  * smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/)
1091  * which improves the precision of the simulation.
1092  *
1093  * @code
1094  * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
1095  * {
1096  * const double layer_transition_center =
1097  * material_a_wavelength / 2 +
1098  * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1099  * material_b_wavelength / 4;
1100  * if (std::abs(p[0]) >=
1101  * (layer_transition_center - average_rho_width / 2) &&
1102  * std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1103  * {
1104  * const double coefficient =
1105  * (std::abs(p[0]) -
1106  * (layer_transition_center - average_rho_width / 2)) /
1107  * average_rho_width;
1108  * return (1 - coefficient) * material_b_rho +
1109  * coefficient * material_a_rho;
1110  * }
1111  * }
1112  *
1113  * @endcode
1114  *
1115  * then the cavity
1116  *
1117  * @code
1118  * if (std::abs(p[0]) <= material_a_wavelength / 2)
1119  * {
1120  * return material_a_rho;
1121  * }
1122  *
1123  * @endcode
1124  *
1125  * the material_a layers
1126  *
1127  * @code
1128  * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
1129  * {
1130  * const double layer_center =
1131  * material_a_wavelength / 2 +
1132  * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1133  * material_b_wavelength / 4 + material_a_wavelength / 8;
1134  * const double layer_width = material_a_wavelength / 4;
1135  * if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1136  * std::abs(p[0]) <= (layer_center + layer_width / 2))
1137  * {
1138  * return material_a_rho;
1139  * }
1140  * }
1141  *
1142  * @endcode
1143  *
1144  * the material_b layers
1145  *
1146  * @code
1147  * for (unsigned int idx = 0; idx < nb_mirror_pairs; idx++)
1148  * {
1149  * const double layer_center =
1150  * material_a_wavelength / 2 +
1151  * idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1152  * material_b_wavelength / 8;
1153  * const double layer_width = material_b_wavelength / 4;
1154  * if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1155  * std::abs(p[0]) <= (layer_center + layer_width / 2))
1156  * {
1157  * return material_b_rho;
1158  * }
1159  * }
1160  *
1161  * @endcode
1162  *
1163  * and finally the default is material_a.
1164  *
1165  * @code
1166  * return material_a_rho;
1167  * }
1168  *
1169  *
1170  *
1171  * @endcode
1172  *
1173  *
1174  * <a name="TheParametersclassimplementation"></a>
1175  * <h4>The `Parameters` class implementation</h4>
1176  *
1177 
1178  *
1179  * The constructor reads all the parameters from the HDF5::Group `data` using
1180  * the HDF5::Group::get_attribute() function.
1181  *
1182  * @code
1183  * template <int dim>
1184  * Parameters<dim>::Parameters(HDF5::Group &data)
1185  * : data(data)
1186  * , simulation_name(data.get_attribute<std::string>("simulation_name"))
1187  * , save_vtu_files(data.get_attribute<bool>("save_vtu_files"))
1188  * , start_frequency(data.get_attribute<double>("start_frequency"))
1189  * , stop_frequency(data.get_attribute<double>("stop_frequency"))
1190  * , nb_frequency_points(data.get_attribute<int>("nb_frequency_points"))
1191  * , lambda(data.get_attribute<double>("lambda"))
1192  * , mu(data.get_attribute<double>("mu"))
1193  * , dimension_x(data.get_attribute<double>("dimension_x"))
1194  * , dimension_y(data.get_attribute<double>("dimension_y"))
1195  * , nb_probe_points(data.get_attribute<int>("nb_probe_points"))
1196  * , grid_level(data.get_attribute<int>("grid_level"))
1197  * , probe_start_point(data.get_attribute<double>("probe_pos_x"),
1198  * data.get_attribute<double>("probe_pos_y") -
1199  * data.get_attribute<double>("probe_width_y") / 2)
1200  * , probe_stop_point(data.get_attribute<double>("probe_pos_x"),
1201  * data.get_attribute<double>("probe_pos_y") +
1202  * data.get_attribute<double>("probe_width_y") / 2)
1203  * , right_hand_side(data)
1204  * , pml(data)
1205  * , rho(data)
1206  * {}
1207  *
1208  *
1209  *
1210  * @endcode
1211  *
1212  *
1213  * <a name="TheQuadratureCacheclassimplementation"></a>
1214  * <h4>The `QuadratureCache` class implementation</h4>
1215  *
1216 
1217  *
1218  * We need to reserve enough space for the mass and stiffness matrices and the
1219  * right hand side vector.
1220  *
1221  * @code
1222  * template <int dim>
1223  * QuadratureCache<dim>::QuadratureCache(const unsigned int dofs_per_cell)
1224  * : dofs_per_cell(dofs_per_cell)
1225  * , mass_coefficient(dofs_per_cell, dofs_per_cell)
1226  * , stiffness_coefficient(dofs_per_cell, dofs_per_cell)
1227  * , right_hand_side(dofs_per_cell)
1228  * {}
1229  *
1230  *
1231  *
1232  * @endcode
1233  *
1234  *
1235  * <a name="ImplementationoftheElasticWaveclass"></a>
1236  * <h3>Implementation of the `ElasticWave` class</h3>
1237  *
1238 
1239  *
1240  *
1241  * <a name="Constructor"></a>
1242  * <h4>Constructor</h4>
1243  *
1244 
1245  *
1246  * This is very similar to the constructor of @ref step_40 "step-40". In addition we create
1247  * the HDF5 datasets `frequency_dataset`, `position_dataset` and
1248  * `displacement`. Note the use of the `template` keyword for the creation of
1249  * the HDF5 datasets. It is a C++ requirement to use the `template` keyword in
1250  * order to treat `create_dataset` as a dependent template name.
1251  *
1252  * @code
1253  * template <int dim>
1254  * ElasticWave<dim>::ElasticWave(const Parameters<dim> &parameters)
1255  * : parameters(parameters)
1256  * , mpi_communicator(MPI_COMM_WORLD)
1257  * , triangulation(mpi_communicator,
1258  * typename Triangulation<dim>::MeshSmoothing(
1259  * Triangulation<dim>::smoothing_on_refinement |
1260  * Triangulation<dim>::smoothing_on_coarsening))
1261  * , quadrature_formula(2)
1262  * , fe(FE_Q<dim>(1), dim)
1263  * , dof_handler(triangulation)
1264  * , frequency(parameters.nb_frequency_points)
1265  * , probe_positions(parameters.nb_probe_points, dim)
1266  * , frequency_dataset(parameters.data.template create_dataset<double>(
1267  * "frequency",
1268  * std::vector<hsize_t>{parameters.nb_frequency_points}))
1269  * , probe_positions_dataset(parameters.data.template create_dataset<double>(
1270  * "position",
1271  * std::vector<hsize_t>{parameters.nb_probe_points, dim}))
1272  * , displacement(
1273  * parameters.data.template create_dataset<std::complex<double>>(
1274  * "displacement",
1275  * std::vector<hsize_t>{parameters.nb_probe_points,
1276  * parameters.nb_frequency_points}))
1277  * , pcout(std::cout,
1278  * (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
1279  * , computing_timer(mpi_communicator,
1280  * pcout,
1281  * TimerOutput::summary,
1282  * TimerOutput::wall_times)
1283  * {}
1284  *
1285  *
1286  *
1287  * @endcode
1288  *
1289  *
1290  * <a name="ElasticWavesetup_system"></a>
1291  * <h4>ElasticWave::setup_system</h4>
1292  *
1293 
1294  *
1295  * There is nothing new in this function, the only difference with @ref step_40 "step-40" is
1296  * that we don't have to apply boundary conditions because we use the PMLs to
1297  * truncate the domain.
1298  *
1299  * @code
1300  * template <int dim>
1301  * void ElasticWave<dim>::setup_system()
1302  * {
1303  * TimerOutput::Scope t(computing_timer, "setup");
1304  *
1305  * dof_handler.distribute_dofs(fe);
1306  *
1307  * locally_owned_dofs = dof_handler.locally_owned_dofs();
1308  * DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
1309  *
1310  * locally_relevant_solution.reinit(locally_owned_dofs,
1311  * locally_relevant_dofs,
1312  * mpi_communicator);
1313  *
1314  * system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1315  *
1316  * constraints.clear();
1317  * constraints.reinit(locally_relevant_dofs);
1318  * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1319  *
1320  * constraints.close();
1321  *
1322  * DynamicSparsityPattern dsp(locally_relevant_dofs);
1323  *
1324  * DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
1326  * locally_owned_dofs,
1327  * mpi_communicator,
1328  * locally_relevant_dofs);
1329  *
1330  * system_matrix.reinit(locally_owned_dofs,
1331  * locally_owned_dofs,
1332  * dsp,
1333  * mpi_communicator);
1334  * }
1335  *
1336  *
1337  *
1338  * @endcode
1339  *
1340  *
1341  * <a name="ElasticWaveassemble_system"></a>
1342  * <h4>ElasticWave::assemble_system</h4>
1343  *
1344 
1345  *
1346  * This function is also very similar to @ref step_40 "step-40", though there are notable
1347  * differences. We assemble the system for each frequency/omega step. In the
1348  * first step we set `calculate_quadrature_data = True` and we calculate the
1349  * mass and stiffness matrices and the right hand side vector. In the
1350  * subsequent steps we will use that data to accelerate the calculation.
1351  *
1352  * @code
1353  * template <int dim>
1354  * void ElasticWave<dim>::assemble_system(const double omega,
1355  * const bool calculate_quadrature_data)
1356  * {
1357  * TimerOutput::Scope t(computing_timer, "assembly");
1358  *
1359  * FEValues<dim> fe_values(fe,
1360  * quadrature_formula,
1363  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
1364  * const unsigned int n_q_points = quadrature_formula.size();
1365  *
1366  * FullMatrix<std::complex<double>> cell_matrix(dofs_per_cell, dofs_per_cell);
1367  * Vector<std::complex<double>> cell_rhs(dofs_per_cell);
1368  *
1369  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1370  *
1371  * @endcode
1372  *
1373  * Here we store the value of the right hand side, rho and the PML.
1374  *
1375  * @code
1376  * std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim));
1377  * std::vector<double> rho_values(n_q_points);
1378  * std::vector<Vector<std::complex<double>>> pml_values(
1379  * n_q_points, Vector<std::complex<double>>(dim));
1380  *
1381  * @endcode
1382  *
1383  * We calculate the stiffness tensor for the @f$\lambda@f$ and @f$\mu@f$ that have
1384  * been defined in the jupyter notebook. Note that contrary to @f$\rho@f$ the
1385  * stiffness is constant among for the whole domain.
1386  *
1387  * @code
1388  * const SymmetricTensor<4, dim> stiffness_tensor =
1389  * get_stiffness_tensor<dim>(parameters.lambda, parameters.mu);
1390  *
1391  * @endcode
1392  *
1393  * We use the same method of @ref step_20 "step-20" for vector-valued problems.
1394  *
1395  * @code
1396  * const FEValuesExtractors::Vector displacement(0);
1397  *
1398  * for (const auto &cell : dof_handler.active_cell_iterators())
1399  * if (cell->is_locally_owned())
1400  * {
1401  * cell_matrix = 0;
1402  * cell_rhs = 0;
1403  *
1404  * @endcode
1405  *
1406  * We have to calculate the values of the right hand side, rho and
1407  * the PML only if we are going to calculate the mass and the
1408  * stiffness matrices. Otherwise we can skip this calculation which
1409  * considerably reduces the total calculation time.
1410  *
1411  * @code
1412  * if (calculate_quadrature_data)
1413  * {
1414  * fe_values.reinit(cell);
1415  *
1416  * parameters.right_hand_side.vector_value_list(
1417  * fe_values.get_quadrature_points(), rhs_values);
1418  * parameters.rho.value_list(fe_values.get_quadrature_points(),
1419  * rho_values);
1420  * parameters.pml.vector_value_list(
1421  * fe_values.get_quadrature_points(), pml_values);
1422  * }
1423  *
1424  * @endcode
1425  *
1426  * We have done this in @ref step_18 "step-18". Get a pointer to the quadrature
1427  * cache data local to the present cell, and, as a defensive
1428  * measure, make sure that this pointer is within the bounds of the
1429  * global array:
1430  *
1431  * @code
1432  * QuadratureCache<dim> *local_quadrature_points_data =
1433  * reinterpret_cast<QuadratureCache<dim> *>(cell->user_pointer());
1434  * Assert(local_quadrature_points_data >= &quadrature_cache.front(),
1435  * ExcInternalError());
1436  * Assert(local_quadrature_points_data <= &quadrature_cache.back(),
1437  * ExcInternalError());
1438  * for (unsigned int q = 0; q < n_q_points; ++q)
1439  * {
1440  * @endcode
1441  *
1442  * The quadrature_data variable is used to store the mass and
1443  * stiffness matrices, the right hand side vector and the value
1444  * of `JxW`.
1445  *
1446  * @code
1447  * QuadratureCache<dim> &quadrature_data =
1448  * local_quadrature_points_data[q];
1449  *
1450  * @endcode
1451  *
1452  * Below we declare the force vector and the parameters of the
1453  * PML @f$s@f$ and @f$\xi@f$.
1454  *
1455  * @code
1456  * Tensor<1, dim> force;
1458  * std::complex<double> xi(1, 0);
1459  *
1460  * @endcode
1461  *
1462  * The following block is calculated only in the first frequency
1463  * step.
1464  *
1465  * @code
1466  * if (calculate_quadrature_data)
1467  * {
1468  * @endcode
1469  *
1470  * Store the value of `JxW`.
1471  *
1472  * @code
1473  * quadrature_data.JxW = fe_values.JxW(q);
1474  *
1475  * for (unsigned int component = 0; component < dim; ++component)
1476  * {
1477  * @endcode
1478  *
1479  * Convert vectors to tensors and calculate xi
1480  *
1481  * @code
1482  * force[component] = rhs_values[q][component];
1483  * s[component] = pml_values[q][component];
1484  * xi *= s[component];
1485  * }
1486  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1487  * {
1488  * const Tensor<1, dim> phi_i =
1489  * fe_values[displacement].value(i, q);
1490  * const Tensor<2, dim> grad_phi_i =
1491  * fe_values[displacement].gradient(i, q);
1492  *
1493  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1494  * {
1495  * const Tensor<1, dim> phi_j =
1496  * fe_values[displacement].value(j, q);
1497  * const Tensor<2, dim> grad_phi_j =
1498  * fe_values[displacement].gradient(j, q);
1499  *
1500  * @endcode
1501  *
1502  * calculate the values of the mass matrix.
1503  *
1504  * @code
1505  * quadrature_data.mass_coefficient[i][j] =
1506  * rho_values[q] * xi * phi_i * phi_j;
1507  *
1508  * @endcode
1509  *
1510  * Loop over the @f$mnkl@f$ indices of the stiffness
1511  * tensor.
1512  *
1513  * @code
1514  * std::complex<double> stiffness_coefficient = 0;
1515  * for (unsigned int m = 0; m < dim; ++m)
1516  * for (unsigned int n = 0; n < dim; ++n)
1517  * for (unsigned int k = 0; k < dim; ++k)
1518  * for (unsigned int l = 0; l < dim; ++l)
1519  * {
1520  * @endcode
1521  *
1522  * Here we calculate the tensors
1523  * @f$\alpha_{mnkl}@f$ and @f$\beta_{mnkl}@f$
1524  *
1525  * @code
1526  * const std::complex<double> alpha =
1527  * xi * stiffness_tensor[m][n][k][l] /
1528  * (2.0 * s[n] * s[k]);
1529  * const std::complex<double> beta =
1530  * xi * stiffness_tensor[m][n][k][l] /
1531  * (2.0 * s[n] * s[l]);
1532  *
1533  * @endcode
1534  *
1535  * Here we calculate the stiffness matrix.
1536  * Note that the stiffness matrix is not
1537  * symmetric because of the PMLs. We use the
1538  * gradient function (see the
1539  * [documentation](https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html))
1540  * which is a <code>Tensor@<2,dim@></code>.
1541  * The matrix @f$G_{ij}@f$ consists of entries
1542  * @f[
1543  * G_{ij}=
1544  * \frac{\partial\phi_i}{\partial x_j}
1545  * =\partial_j \phi_i
1546  * @f]
1547  * Note the position of the indices @f$i@f$ and
1548  * @f$j@f$ and the notation that we use in this
1549  * tutorial: @f$\partial_j\phi_i@f$. As the
1550  * stiffness tensor is not symmetric, it is
1551  * very easy to make a mistake.
1552  *
1553  * @code
1554  * stiffness_coefficient +=
1555  * grad_phi_i[m][n] *
1556  * (alpha * grad_phi_j[l][k] +
1557  * beta * grad_phi_j[k][l]);
1558  * }
1559  *
1560  * @endcode
1561  *
1562  * We save the value of the stiffness matrix in
1563  * quadrature_data
1564  *
1565  * @code
1566  * quadrature_data.stiffness_coefficient[i][j] =
1567  * stiffness_coefficient;
1568  * }
1569  *
1570  * @endcode
1571  *
1572  * and the value of the right hand side in
1573  * quadrature_data.
1574  *
1575  * @code
1576  * quadrature_data.right_hand_side[i] =
1577  * phi_i * force * fe_values.JxW(q);
1578  * }
1579  * }
1580  *
1581  * @endcode
1582  *
1583  * We loop again over the degrees of freedom of the cells to
1584  * calculate the system matrix. These loops are really quick
1585  * because we have already calculated the stiffness and mass
1586  * matrices, only the value of @f$\omega@f$ changes.
1587  *
1588  * @code
1589  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1590  * {
1591  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1592  * {
1593  * std::complex<double> matrix_sum = 0;
1594  * matrix_sum += -std::pow(omega, 2) *
1595  * quadrature_data.mass_coefficient[i][j];
1596  * matrix_sum += quadrature_data.stiffness_coefficient[i][j];
1597  * cell_matrix(i, j) += matrix_sum * quadrature_data.JxW;
1598  * }
1599  * cell_rhs(i) += quadrature_data.right_hand_side[i];
1600  * }
1601  * }
1602  * cell->get_dof_indices(local_dof_indices);
1603  * constraints.distribute_local_to_global(cell_matrix,
1604  * cell_rhs,
1605  * local_dof_indices,
1606  * system_matrix,
1607  * system_rhs);
1608  * }
1609  *
1610  * system_matrix.compress(VectorOperation::add);
1611  * system_rhs.compress(VectorOperation::add);
1612  * }
1613  *
1614  * @endcode
1615  *
1616  *
1617  * <a name="ElasticWavesolve"></a>
1618  * <h4>ElasticWave::solve</h4>
1619  *
1620 
1621  *
1622  * This is even more simple than in @ref step_40 "step-40". We use the parallel direct solver
1623  * MUMPS which requires less options than an iterative solver. The drawback is
1624  * that it does not scale very well. It is not straightforward to solve the
1625  * Helmholtz equation with an iterative solver. The shifted Laplacian
1626  * multigrid method is a well known approach to precondition this system, but
1627  * this is beyond the scope of this tutorial.
1628  *
1629  * @code
1630  * template <int dim>
1631  * void ElasticWave<dim>::solve()
1632  * {
1633  * TimerOutput::Scope t(computing_timer, "solve");
1634  * LinearAlgebraPETSc::MPI::Vector completely_distributed_solution(
1635  * locally_owned_dofs, mpi_communicator);
1636  *
1637  * SolverControl solver_control;
1638  * PETScWrappers::SparseDirectMUMPS solver(solver_control, mpi_communicator);
1639  * solver.solve(system_matrix, completely_distributed_solution, system_rhs);
1640  *
1641  * pcout << " Solved in " << solver_control.last_step() << " iterations."
1642  * << std::endl;
1643  * constraints.distribute(completely_distributed_solution);
1644  * locally_relevant_solution = completely_distributed_solution;
1645  * }
1646  *
1647  * @endcode
1648  *
1649  *
1650  * <a name="ElasticWaveinitialize_position_vector"></a>
1651  * <h4>ElasticWave::initialize_position_vector</h4>
1652  *
1653 
1654  *
1655  * We use this function to calculate the values of the position vector.
1656  *
1657  * @code
1658  * template <int dim>
1659  * void ElasticWave<dim>::initialize_probe_positions_vector()
1660  * {
1661  * for (unsigned int position_idx = 0;
1662  * position_idx < parameters.nb_probe_points;
1663  * ++position_idx)
1664  * {
1665  * @endcode
1666  *
1667  * Because of the way the operator + and - are overloaded to subtract
1668  * two points, the following has to be done:
1669  * `Point_b<dim> + (-Point_a<dim>)`
1670  *
1671  * @code
1672  * const Point<dim> p =
1673  * (position_idx / ((double)(parameters.nb_probe_points - 1))) *
1674  * (parameters.probe_stop_point + (-parameters.probe_start_point)) +
1675  * parameters.probe_start_point;
1676  * probe_positions[position_idx][0] = p[0];
1677  * probe_positions[position_idx][1] = p[1];
1678  * if (dim == 3)
1679  * {
1680  * probe_positions[position_idx][2] = p[2];
1681  * }
1682  * }
1683  * }
1684  *
1685  * @endcode
1686  *
1687  *
1688  * <a name="ElasticWavestore_frequency_step_data"></a>
1689  * <h4>ElasticWave::store_frequency_step_data</h4>
1690  *
1691 
1692  *
1693  * This function stores in the HDF5 file the measured energy by the probe.
1694  *
1695  * @code
1696  * template <int dim>
1697  * void
1698  * ElasticWave<dim>::store_frequency_step_data(const unsigned int frequency_idx)
1699  * {
1700  * TimerOutput::Scope t(computing_timer, "store_frequency_step_data");
1701  *
1702  * @endcode
1703  *
1704  * We store the displacement in the @f$x@f$ direction; the displacement in the
1705  * @f$y@f$ direction is negligible.
1706  *
1707  * @code
1708  * const unsigned int probe_displacement_component = 0;
1709  *
1710  * @endcode
1711  *
1712  * The vector coordinates contains the coordinates in the HDF5 file of the
1713  * points of the probe that are located in locally owned cells. The vector
1714  * displacement_data contains the value of the displacement at these points.
1715  *
1716  * @code
1717  * std::vector<hsize_t> coordinates;
1718  * std::vector<std::complex<double>> displacement_data;
1719  * for (unsigned int position_idx = 0;
1720  * position_idx < parameters.nb_probe_points;
1721  * ++position_idx)
1722  * {
1723  * Point<dim> point;
1724  * for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1725  * {
1726  * point[dim_idx] = probe_positions[position_idx][dim_idx];
1727  * }
1728  * bool point_in_locally_owned_cell;
1729  * {
1730  * @endcode
1731  *
1732  * First we have to find out if the point is in a locally owned cell.
1733  *
1734  * @code
1735  * auto mapping = StaticMappingQ1<dim>::mapping;
1736  * const std::pair<typename DoFHandler<dim>::active_cell_iterator,
1737  * Point<dim>>
1738  * cell_point = GridTools::find_active_cell_around_point(mapping,
1739  * dof_handler,
1740  * point);
1741  *
1742  * point_in_locally_owned_cell = cell_point.first->is_locally_owned();
1743  * }
1744  * if (point_in_locally_owned_cell)
1745  * {
1746  * @endcode
1747  *
1748  * Then we can store the values of the displacement in the points of
1749  * the probe in `displacement_data`.
1750  *
1751  * @code
1752  * Vector<std::complex<double>> tmp_vector(dim);
1753  * VectorTools::point_value(dof_handler,
1754  * locally_relevant_solution,
1755  * point,
1756  * tmp_vector);
1757  * coordinates.emplace_back(position_idx);
1758  * coordinates.emplace_back(frequency_idx);
1759  * displacement_data.emplace_back(
1760  * tmp_vector(probe_displacement_component));
1761  * }
1762  * }
1763  *
1764  * @endcode
1765  *
1766  * We write the displacement data in the HDF5 file. The call
1767  * HDF5::DataSet::write_selection() is MPI collective which means that all
1768  * the processes have to participate.
1769  *
1770  * @code
1771  * if (coordinates.size() > 0)
1772  * {
1773  * displacement.write_selection(displacement_data, coordinates);
1774  * }
1775  * @endcode
1776  *
1777  * Therefore even if the process has no data to write it has to participate
1778  * in the collective call. For this we can use HDF5::DataSet::write_none().
1779  * Note that we have to specify the data type, in this case
1780  * `std::complex<double>`.
1781  *
1782  * @code
1783  * else
1784  * {
1785  * displacement.write_none<std::complex<double>>();
1786  * }
1787  *
1788  * @endcode
1789  *
1790  * If the variable `save_vtu_files` in the input file equals `True` then all
1791  * the data will be saved as vtu. The procedure to write `vtu` files has
1792  * been described in @ref step_40 "step-40".
1793  *
1794  * @code
1795  * if (parameters.save_vtu_files)
1796  * {
1797  * std::vector<std::string> solution_names(dim, "displacement");
1798  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1799  * interpretation(
1801  *
1802  * DataOut<dim> data_out;
1803  * data_out.add_data_vector(dof_handler,
1804  * locally_relevant_solution,
1805  * solution_names,
1806  * interpretation);
1807  * Vector<float> subdomain(triangulation.n_active_cells());
1808  * for (unsigned int i = 0; i < subdomain.size(); ++i)
1809  * subdomain(i) = triangulation.locally_owned_subdomain();
1810  * data_out.add_data_vector(subdomain, "subdomain");
1811  *
1812  * std::vector<Vector<double>> force(
1813  * dim, Vector<double>(triangulation.n_active_cells()));
1814  * std::vector<Vector<double>> pml(
1815  * dim, Vector<double>(triangulation.n_active_cells()));
1816  * Vector<double> rho(triangulation.n_active_cells());
1817  *
1818  * for (auto &cell : triangulation.active_cell_iterators())
1819  * {
1820  * if (cell->is_locally_owned())
1821  * {
1822  * for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1823  * {
1824  * force[dim_idx](cell->active_cell_index()) =
1825  * parameters.right_hand_side.value(cell->center(), dim_idx);
1826  * pml[dim_idx](cell->active_cell_index()) =
1827  * parameters.pml.value(cell->center(), dim_idx).imag();
1828  * }
1829  * rho(cell->active_cell_index()) =
1830  * parameters.rho.value(cell->center());
1831  * }
1832  * @endcode
1833  *
1834  * And on the cells that we are not interested in, set the
1835  * respective value to a bogus value in order to make sure that if
1836  * we were somehow wrong about our assumption we would find out by
1837  * looking at the graphical output:
1838  *
1839  * @code
1840  * else
1841  * {
1842  * for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1843  * {
1844  * force[dim_idx](cell->active_cell_index()) = -1e+20;
1845  * pml[dim_idx](cell->active_cell_index()) = -1e+20;
1846  * }
1847  * rho(cell->active_cell_index()) = -1e+20;
1848  * }
1849  * }
1850  *
1851  * for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1852  * {
1853  * data_out.add_data_vector(force[dim_idx],
1854  * "force_" + std::to_string(dim_idx));
1855  * data_out.add_data_vector(pml[dim_idx],
1856  * "pml_" + std::to_string(dim_idx));
1857  * }
1858  * data_out.add_data_vector(rho, "rho");
1859  *
1860  * data_out.build_patches();
1861  *
1862  * std::stringstream frequency_idx_stream;
1863  * const unsigned int nb_number_positions =
1864  * ((unsigned int)std::log10(parameters.nb_frequency_points)) + 1;
1865  * frequency_idx_stream << std::setw(nb_number_positions)
1866  * << std::setfill('0') << frequency_idx;
1867  * std::string filename = (parameters.simulation_name + "_" +
1868  * frequency_idx_stream.str() + ".vtu");
1869  * data_out.write_vtu_in_parallel(filename.c_str(), mpi_communicator);
1870  * }
1871  * }
1872  *
1873  *
1874  *
1875  * @endcode
1876  *
1877  *
1878  * <a name="ElasticWaveoutput_results"></a>
1879  * <h4>ElasticWave::output_results</h4>
1880  *
1881 
1882  *
1883  * This function writes the datasets that have not already been written.
1884  *
1885  * @code
1886  * template <int dim>
1887  * void ElasticWave<dim>::output_results()
1888  * {
1889  * @endcode
1890  *
1891  * The vectors `frequency` and `position` are the same for all the
1892  * processes. Therefore any of the processes can write the corresponding
1893  * `datasets`. Because the call HDF5::DataSet::write is MPI collective, the
1894  * rest of the processes will have to call HDF5::DataSet::write_none.
1895  *
1896  * @code
1897  * if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
1898  * {
1899  * frequency_dataset.write(frequency);
1900  * probe_positions_dataset.write(probe_positions);
1901  * }
1902  * else
1903  * {
1904  * frequency_dataset.write_none<double>();
1905  * probe_positions_dataset.write_none<double>();
1906  * }
1907  * }
1908  *
1909  *
1910  *
1911  * @endcode
1912  *
1913  *
1914  * <a name="ElasticWavesetup_quadrature_cache"></a>
1915  * <h4>ElasticWave::setup_quadrature_cache</h4>
1916  *
1917 
1918  *
1919  * We use this function at the beginning of our computations to set up initial
1920  * values of the cache variables. This function has been described in @ref step_18 "step-18".
1921  * There are no differences with the function of @ref step_18 "step-18".
1922  *
1923  * @code
1924  * template <int dim>
1925  * void ElasticWave<dim>::setup_quadrature_cache()
1926  * {
1927  * triangulation.clear_user_data();
1928  *
1929  * {
1930  * std::vector<QuadratureCache<dim>> tmp;
1931  * quadrature_cache.swap(tmp);
1932  * }
1933  *
1934  * quadrature_cache.resize(triangulation.n_locally_owned_active_cells() *
1935  * quadrature_formula.size(),
1936  * QuadratureCache<dim>(fe.dofs_per_cell));
1937  * unsigned int cache_index = 0;
1938  * for (const auto &cell : triangulation.active_cell_iterators())
1939  * if (cell->is_locally_owned())
1940  * {
1941  * cell->set_user_pointer(&quadrature_cache[cache_index]);
1942  * cache_index += quadrature_formula.size();
1943  * }
1944  * Assert(cache_index == quadrature_cache.size(), ExcInternalError());
1945  * }
1946  *
1947  *
1948  *
1949  * @endcode
1950  *
1951  *
1952  * <a name="ElasticWavefrequency_sweep"></a>
1953  * <h4>ElasticWave::frequency_sweep</h4>
1954  *
1955 
1956  *
1957  * For clarity we divide the function `run` of @ref step_40 "step-40" into the functions
1958  * `run` and `frequency_sweep`. In the function `frequency_sweep` we place the
1959  * iteration over the frequency vector.
1960  *
1961  * @code
1962  * template <int dim>
1963  * void ElasticWave<dim>::frequency_sweep()
1964  * {
1965  * for (unsigned int frequency_idx = 0;
1966  * frequency_idx < parameters.nb_frequency_points;
1967  * ++frequency_idx)
1968  * {
1969  * pcout << parameters.simulation_name + " frequency idx: "
1970  * << frequency_idx << '/' << parameters.nb_frequency_points - 1
1971  * << std::endl;
1972  *
1973  *
1974  *
1975  * setup_system();
1976  * if (frequency_idx == 0)
1977  * {
1978  * pcout << " Number of active cells : "
1979  * << triangulation.n_active_cells() << std::endl;
1980  * pcout << " Number of degrees of freedom : "
1981  * << dof_handler.n_dofs() << std::endl;
1982  * }
1983  *
1984  * if (frequency_idx == 0)
1985  * {
1986  * @endcode
1987  *
1988  * Write the simulation parameters only once
1989  *
1990  * @code
1991  * parameters.data.set_attribute("active_cells",
1992  * triangulation.n_active_cells());
1993  * parameters.data.set_attribute("degrees_of_freedom",
1994  * dof_handler.n_dofs());
1995  * }
1996  *
1997  * @endcode
1998  *
1999  * We calculate the frequency and omega values for this particular step.
2000  *
2001  * @code
2002  * const double current_loop_frequency =
2003  * (parameters.start_frequency +
2004  * frequency_idx *
2005  * (parameters.stop_frequency - parameters.start_frequency) /
2006  * (parameters.nb_frequency_points - 1));
2007  * const double current_loop_omega =
2008  * 2 * numbers::PI * current_loop_frequency;
2009  *
2010  * @endcode
2011  *
2012  * In the first frequency step we calculate the mass and stiffness
2013  * matrices and the right hand side. In the subsequent frequency steps
2014  * we will use those values. This improves considerably the calculation
2015  * time.
2016  *
2017  * @code
2018  * assemble_system(current_loop_omega,
2019  * (frequency_idx == 0) ? true : false);
2020  * solve();
2021  *
2022  * frequency[frequency_idx] = current_loop_frequency;
2023  * store_frequency_step_data(frequency_idx);
2024  *
2025  * computing_timer.print_summary();
2026  * computing_timer.reset();
2027  * pcout << std::endl;
2028  * }
2029  * }
2030  *
2031  *
2032  *
2033  * @endcode
2034  *
2035  *
2036  * <a name="ElasticWaverun"></a>
2037  * <h4>ElasticWave::run</h4>
2038  *
2039 
2040  *
2041  * This function is very similar to the one in @ref step_40 "step-40".
2042  *
2043  * @code
2044  * template <int dim>
2045  * void ElasticWave<dim>::run()
2046  * {
2047  * #ifdef DEBUG
2048  * pcout << "Debug mode" << std::endl;
2049  * #else
2050  * pcout << "Release mode" << std::endl;
2051  * #endif
2052  *
2053  * {
2054  * Point<dim> p1;
2055  * p1(0) = -parameters.dimension_x / 2;
2056  * p1(1) = -parameters.dimension_y / 2;
2057  * if (dim == 3)
2058  * {
2059  * p1(2) = -parameters.dimension_y / 2;
2060  * }
2061  * Point<dim> p2;
2062  * p2(0) = parameters.dimension_x / 2;
2063  * p2(1) = parameters.dimension_y / 2;
2064  * if (dim == 3)
2065  * {
2066  * p2(2) = parameters.dimension_y / 2;
2067  * }
2068  * std::vector<unsigned int> divisions(dim);
2069  * divisions[0] = int(parameters.dimension_x / parameters.dimension_y);
2070  * divisions[1] = 1;
2071  * if (dim == 3)
2072  * {
2073  * divisions[2] = 1;
2074  * }
2076  * divisions,
2077  * p1,
2078  * p2);
2079  * }
2080  *
2081  * triangulation.refine_global(parameters.grid_level);
2082  *
2083  * setup_quadrature_cache();
2084  *
2085  * initialize_probe_positions_vector();
2086  *
2087  * frequency_sweep();
2088  *
2089  * output_results();
2090  * }
2091  * } // namespace step62
2092  *
2093  *
2094  *
2095  * @endcode
2096  *
2097  *
2098  * <a name="Themainfunction"></a>
2099  * <h4>The main function</h4>
2100  *
2101 
2102  *
2103  * The main function is very similar to the one in @ref step_40 "step-40".
2104  *
2105  * @code
2106  * int main(int argc, char *argv[])
2107  * {
2108  * try
2109  * {
2110  * using namespace dealii;
2111  * const unsigned int dim = 2;
2112  *
2113  * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
2114  *
2115  * HDF5::File data_file("results.h5",
2117  * MPI_COMM_WORLD);
2118  * HDF5::Group data = data_file.open_group("data");
2119  *
2120  * {
2121  * @endcode
2122  *
2123  * Displacement simulation. The parameters are read from the
2124  * displacement HDF5 group and the results are saved in the same HDF5
2125  * group.
2126  *
2127  * @code
2128  * auto displacement = data.open_group("displacement");
2129  * step62::Parameters<dim> parameters(displacement);
2130  *
2131  * step62::ElasticWave<dim> elastic_problem(parameters);
2132  * elastic_problem.run();
2133  * }
2134  *
2135  * {
2136  * @endcode
2137  *
2138  * Calibration simulation. The parameters are read from the calibration
2139  * HDF5 group and the results are saved in the same HDF5 group.
2140  *
2141  * @code
2142  * auto calibration = data.open_group("calibration");
2143  * step62::Parameters<dim> parameters(calibration);
2144  *
2145  * step62::ElasticWave<dim> elastic_problem(parameters);
2146  * elastic_problem.run();
2147  * }
2148  * }
2149  * catch (std::exception &exc)
2150  * {
2151  * std::cerr << std::endl
2152  * << std::endl
2153  * << "----------------------------------------------------"
2154  * << std::endl;
2155  * std::cerr << "Exception on processing: " << std::endl
2156  * << exc.what() << std::endl
2157  * << "Aborting!" << std::endl
2158  * << "----------------------------------------------------"
2159  * << std::endl;
2160  *
2161  * return 1;
2162  * }
2163  * catch (...)
2164  * {
2165  * std::cerr << std::endl
2166  * << std::endl
2167  * << "----------------------------------------------------"
2168  * << std::endl;
2169  * std::cerr << "Unknown exception!" << std::endl
2170  * << "Aborting!" << std::endl
2171  * << "----------------------------------------------------"
2172  * << std::endl;
2173  * return 1;
2174  * }
2175  *
2176  * return 0;
2177  * }
2178  * @endcode
2179 <a name="Results"></a><h1>Results</h1>
2180 
2181 
2182 <a name="Resonancefrequencyandbandgap"></a><h3>Resonance frequency and bandgap</h3>
2183 
2184 
2185 The results are analyzed in the
2186 [jupyter notebook](https://github.com/dealii/dealii/blob/phononic-cavity/examples/step-62/step-62.ipynb)
2187 with the following code
2188 @code{.py}
2189 h5_file = h5py.File('results.h5', 'r')
2190 data = h5_file['data']
2191 
2192 # Gaussian function that we use to fit the resonance
2193 def resonance_f(freq, freq_m, quality_factor, max_amplitude):
2194  omega = 2 * constants.pi * freq
2195  omega_m = 2 * constants.pi * freq_m
2196  gamma = omega_m / quality_factor
2197  return max_amplitude * omega_m**2 * gamma**2 / (((omega_m**2 - omega**2)**2 + gamma**2 * omega**2))
2198 
2199 frequency = data['displacement']['frequency'][...]
2200 # Average the probe points
2201 displacement = np.mean(data['displacement']['displacement'], axis=0)
2202 calibration_displacement = np.mean(data['calibration']['displacement'], axis=0)
2203 reflection_coefficient = displacement / calibration_displacement
2204 reflectivity = (np.abs(np.mean(data['displacement']['displacement'][...]**2, axis=0))/
2205  np.abs(np.mean(data['calibration']['displacement'][...]**2, axis=0)))
2206 
2207 try:
2208  x_data = frequency
2209  y_data = reflectivity
2210  quality_factor_guess = 1e3
2211  freq_guess = x_data[np.argmax(y_data)]
2212  amplitude_guess = np.max(y_data)
2213  fit_result, covariance = scipy.optimize.curve_fit(resonance_f, x_data, y_data,
2214  [freq_guess, quality_factor_guess, amplitude_guess])
2215  freq_m = fit_result[0]
2216  quality_factor = np.abs(fit_result[1])
2217  max_amplitude = fit_result[2]
2218  y_data_fit = resonance_f(x_data, freq_m, quality_factor, max_amplitude)
2219 
2220  fig = plt.figure()
2221  plt.plot(frequency / 1e9, reflectivity, frequency / 1e9, y_data_fit)
2222  plt.xlabel('frequency (GHz)')
2223  plt.ylabel('amplitude^2 (a.u.)')
2224  plt.title('Transmission\n' + 'freq = ' + "%.7g" % (freq_guess / 1e9) + 'GHz Q = ' + "%.6g" % quality_factor)
2225 except:
2226  fig = plt.figure()
2227  plt.plot(frequency / 1e9, reflectivity)
2228  plt.xlabel('frequency (GHz)')
2229  plt.ylabel('amplitude^2 (a.u.)')
2230  plt.title('Transmission')
2231 
2232 fig = plt.figure()
2233 plt.plot(frequency / 1e9, np.angle(reflection_coefficient))
2234 plt.xlabel('frequency (GHz)')
2235 plt.ylabel('phase (rad)')
2236 plt.title('Phase (transmission coefficient)\n')
2237 
2238 plt.show()
2239 h5_file.close()
2240 @endcode
2241 
2242 A phononic cavity is characterized by the
2243 [resonance frequency](https://en.wikipedia.org/wiki/Resonance) and the
2244 [the quality factor](https://en.wikipedia.org/wiki/Q_factor).
2245 The quality factor is equal to the ratio between the stored energy in the resonator and the energy
2246 dissipated energy per cycle, which is approximately equivalent to the ratio between the
2247 resonance frequency and the
2248 [full width at half maximum (FWHM)](https://en.wikipedia.org/wiki/Full_width_at_half_maximum).
2249 The FWHM is equal to the bandwidth over which the power of vibration is greater than half the
2250 power at the resonant frequency.
2251 @f[
2252 Q = \frac{f_r}{\Delta f} = \frac{\omega_r}{\Delta \omega} =
2253 2 \pi \times \frac{\text{energy stored}}{\text{energy dissipated per cycle}}
2254 @f]
2255 
2256 The square of the amplitude of the mechanical resonance @f$a^2@f$ as a function of the frequency
2257 has a gaussian shape
2258 @f[
2259 a^2 = a_\textrm{max}^2\frac{\omega^2\Gamma^2}{(\omega_r^2-\omega^2)^2+\Gamma^2\omega^2}
2260 @f]
2261 where @f$f_r = \frac{\omega_r}{2\pi}@f$ is the resonance frequency and @f$\Gamma=\frac{\omega_r}{Q}@f$ is the dissipation rate.
2262 We used the previous equation in the jupyter notebook to fit the mechanical resonance.
2263 
2264 Given the values we have chosen for the parameters, one could estimate the resonance frequency
2265 analytically. Indeed, this is then confirmed by what we get in this program:
2266 the phononic superlattice cavity exhibits a mechanical resonance at 20GHz and a quality factor of 5046.
2267 The following images show the transmission amplitude and phase as a function of frequency in the
2268 vicinity of the resonance frequency:
2269 
2270 <img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.05.png" height="400" />
2271 <img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.06.png" height="400" />
2272 
2273 The images above suggest that the periodic structure has its intended effect: It really only lets waves of a very
2274 specific frequency pass through, whereas all other waves are reflected. This is of course precisely what one builds
2275 these sorts of devices for.
2276 But it is not quite this easy. In practice, there is really only a "band gap", i.e., the device blocks waves other than
2277 the desired one at 20GHz only within a certain frequency range. Indeed, to find out how large this "gap" is within
2278 which waves are blocked, we can extend the frequency range to 16 GHz through the appropriate parameters in the
2279 input file. We then obtain the following image:
2280 
2281 <img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.07.png" height="400" />
2282 
2283 What this image suggests is that in the range of around 18 to around 22 GHz, really only the waves with a frequency
2284 of 20 GHz are allowed to pass through, but beyond this range, there are plenty of other frequencies that can pass
2285 through the device.
2286 
2287 <a name="Modeprofile"></a><h3>Mode profile</h3>
2288 
2289 
2290 We can inspect the mode profile with Paraview or VisIt.
2291 As we have discussed, at resonance all the mechanical
2292 energy is transmitted and the amplitude of motion is amplified inside the cavity.
2293 It can be observed that the PMLs are quite effective to truncate the solution.
2294 The following image shows the mode profile at resonance:
2295 
2296 <img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.08.png" height="400" />
2297 
2298 On the other hand, out of resonance all the mechanical energy is
2299 reflected. The following image shows the profile at 19.75 GHz.
2300 Note the interference between the force pulse and the reflected wave
2301 at the position @f$x=-8\mu\textrm{m}@f$.
2302 
2303 <img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.09.png" height="400" />
2304 
2305 <a name="Experimentalapplications"></a><h3>Experimental applications</h3>
2306 
2307 
2308 Phononic superlattice cavities find application in
2309 [quantum optomechanics](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.86.1391).
2310 Here we have presented the simulation of a 2D superlattice cavity,
2311 but this code can be used as well to simulate "real world" 3D devices such as
2312 [micropillar superlattice cavities](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.060101),
2313 which are promising candidates to study macroscopic quantum phenomena.
2314 The 20GHz mode of a micropillar superlattice cavity is essentially a mechanical harmonic oscillator that is very well isolated
2315 from the environment. If the device is cooled down to 20mK in a dilution fridge, the mode would then become a
2316 macroscopic quantum harmonic oscillator.
2317  *
2318  *
2319 <a name="PlainProg"></a>
2320 <h1> The plain program</h1>
2321 @include "step-62.cc"
2322 */
HDF5::File::FileAccessMode::open
@ open
internal::QGaussLobatto::gamma
long double gamma(const unsigned int n)
Definition: quadrature_lib.cc:96
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
HDF5::File
Definition: hdf5.h:1067
TimerOutput::Scope
Definition: timer.h:554
HDF5::Group::open_group
Group open_group(const std::string &name) const
Definition: hdf5.cc:1345
StaticMappingQ1
Definition: mapping_q1.h:88
SymmetricTensor< 4, dim >
dealii
Definition: namespace_dealii.h:25
HDF5
Definition: hdf5.h:331
DataOutBase::vtu
@ vtu
Definition: data_out_base.h:1610
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2360
VectorTools::point_value
void point_value(const DoFHandler< dim, spacedim > &dof, const VectorType &fe_function, const Point< spacedim, double > &point, Vector< typename VectorType::value_type > &value)
VectorOperation::add
@ add
Definition: vector_operation.h:53
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
LocalIntegrators::Advection::cell_matrix
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double >> &velocity, const double factor=1.)
Definition: advection.h:80
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
SolverControl::last_step
unsigned int last_step() const
Definition: solver_control.cc:127
FEValues< dim >
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
GridTools::find_active_cell_around_point
MeshType< dim, spacedim >::active_cell_iterator find_active_cell_around_point(const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={})
Definition: grid_tools_dof_handlers.cc:549
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandlerType &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition: dof_tools_sparsity.cc:63
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:115
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
FEValuesExtractors::Vector
Definition: fe_values_extractors.h:150
angle
const double angle
Definition: grid_tools_nontemplates.cc:277
HDF5::DataSet::write
void write(const Container &data)
Definition: hdf5.cc:942
Tensor< 1, dim >
GridGenerator::subdivided_hyper_rectangle
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
VectorTools::mean
@ mean
Definition: vector_tools_common.h:79
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
DynamicSparsityPattern
Definition: dynamic_sparsity_pattern.h:323
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandlerType &dof_handler, AffineConstraints< number > &constraints)
Definition: dof_tools_constraints.cc:1725
PETScWrappers::MPI::Vector
Definition: petsc_vector.h:158
PETScWrappers::SparseDirectMUMPS
Definition: petsc_solver.h:937
SIMDComparison::equal
@ equal
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5450
value
static const bool value
Definition: dof_tools_constraints.cc:433
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
DataOutInterface::write_vtu_in_parallel
void write_vtu_in_parallel(const std::string &filename, MPI_Comm comm) const
Definition: data_out_base.cc:6886
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
DoFTools::extract_locally_relevant_dofs
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1173
std::pow
inline ::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &x, const Number p)
Definition: vectorization.h:5428
SparsityTools::distribute_sparsity_pattern
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm &mpi_comm, const IndexSet &locally_relevant_rows)
Definition: sparsity_tools.cc:1046
HDF5::Group
Definition: hdf5.h:968
Point< dim >
FullMatrix
Definition: full_matrix.h:71
HDF5::DataSet::write_none
void write_none()
Definition: hdf5.cc:1132
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
HDF5::DataSet::write_selection
void write_selection(const Container &data, const std::vector< hsize_t > &coordinates)
Definition: hdf5.cc:974
SolverControl
Definition: solver_control.h:67
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
Differentiation::SD::log10
Expression log10(const Expression &x)
Definition: symengine_math.cc:84
first
Point< 2 > first
Definition: grid_out.cc:4352
DataOut< dim >
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector
Definition: mapping_q1_eulerian.h:32
parallel
Definition: distributed.h:416
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
DataComponentInterpretation::component_is_part_of_vector
@ component_is_part_of_vector
Definition: data_component_interpretation.h:61
int
Utilities::MPI::MPI_InitFinalize
Definition: mpi.h:828
DataOut_DoFData::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
Definition: data_out_dof_data.h:1090