deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
step-62.h
Go to the documentation of this file.
1) const
1019 *   {
1020 * @endcode
1021 *
1023 * @f[
1024 * c = \frac{K_e}{\rho}
1025 * @f]
1028 * than the wavelength. In this case it can be shown that for the two
1029 * dimensional case
1030 * @f[
1031 * K_e = 4\mu\frac{\lambda +\mu}{\lambda+2\mu}
1032 * @f]
1033 * and for the three dimensional case @f$K_e@f$ is equal to the Young's modulus.
1034 * @f[
1035 * K_e = \mu\frac{3\lambda +2\mu}{\lambda+\mu}
1036 * @f]
1037 *
1038 * @code
1039 *   double elastic_constant;
1040 *   if (dim == 2)
1041 *   {
1042 *   elastic_constant = 4 * mu * (lambda + mu) / (lambda + 2 * mu);
1043 *   }
1044 *   else if (dim == 3)
1045 *   {
1046 *   elastic_constant = mu * (3 * lambda + 2 * mu) / (lambda + mu);
1047 *   }
1048 *   else
1049 *   DEAL_II_NOT_IMPLEMENTED();
1050 *  
1051 *   const double material_a_speed_of_sound =
1052 *   std::sqrt(elastic_constant / material_a_rho);
1053 *   const double material_a_wavelength =
1054 *   material_a_speed_of_sound / cavity_resonance_frequency;
1055 *   const double material_b_speed_of_sound =
1056 *   std::sqrt(elastic_constant / material_b_rho);
1057 *   const double material_b_wavelength =
1058 *   material_b_speed_of_sound / cavity_resonance_frequency;
1059 *  
1060 * @endcode
1061 *
1062 * The density @f$\rho@f$ takes the following form
1063 * <img alt="Phononic superlattice cavity"
1064 * src="https://www.dealii.org/images/steps/developer/step-62.04.svg"
1065 * height="200" />
1066 * where the brown color represents material_a and the green color
1067 * represents material_b.
1068 *
1069 * @code
1070 *   for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1071 *   {
1072 *   const double layer_transition_center =
1073 *   material_a_wavelength / 2 +
1074 *   idx * (material_b_wavelength / 4 + material_a_wavelength / 4);
1075 *   if (std::abs(p[0]) >=
1076 *   (layer_transition_center - average_rho_width / 2) &&
1077 *   std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1078 *   {
1079 *   const double coefficient =
1080 *   (std::abs(p[0]) -
1081 *   (layer_transition_center - average_rho_width / 2)) /
1082 *   average_rho_width;
1083 *   return (1 - coefficient) * material_a_rho +
1084 *   coefficient * material_b_rho;
1085 *   }
1086 *   }
1087 *  
1088 * @endcode
1089 *
1090 * Here we define the
1091 * [subpixel
1092 * smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/)
1093 * which improves the precision of the simulation.
1094 *
1095 * @code
1096 *   for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1097 *   {
1098 *   const double layer_transition_center =
1099 *   material_a_wavelength / 2 +
1100 *   idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1101 *   material_b_wavelength / 4;
1102 *   if (std::abs(p[0]) >=
1103 *   (layer_transition_center - average_rho_width / 2) &&
1104 *   std::abs(p[0]) <= (layer_transition_center + average_rho_width / 2))
1105 *   {
1106 *   const double coefficient =
1107 *   (std::abs(p[0]) -
1108 *   (layer_transition_center - average_rho_width / 2)) /
1109 *   average_rho_width;
1110 *   return (1 - coefficient) * material_b_rho +
1111 *   coefficient * material_a_rho;
1112 *   }
1113 *   }
1114 *  
1115 * @endcode
1116 *
1117 * then the cavity
1118 *
1119 * @code
1120 *   if (std::abs(p[0]) <= material_a_wavelength / 2)
1121 *   {
1122 *   return material_a_rho;
1123 *   }
1124 *  
1125 * @endcode
1126 *
1127 * the material_a layers
1128 *
1129 * @code
1130 *   for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1131 *   {
1132 *   const double layer_center =
1133 *   material_a_wavelength / 2 +
1134 *   idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1135 *   material_b_wavelength / 4 + material_a_wavelength / 8;
1136 *   const double layer_width = material_a_wavelength / 4;
1137 *   if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1138 *   std::abs(p[0]) <= (layer_center + layer_width / 2))
1139 *   {
1140 *   return material_a_rho;
1141 *   }
1142 *   }
1143 *  
1144 * @endcode
1145 *
1146 * the material_b layers
1147 *
1148 * @code
1149 *   for (unsigned int idx = 0; idx < nb_mirror_pairs; ++idx)
1150 *   {
1151 *   const double layer_center =
1152 *   material_a_wavelength / 2 +
1153 *   idx * (material_b_wavelength / 4 + material_a_wavelength / 4) +
1154 *   material_b_wavelength / 8;
1155 *   const double layer_width = material_b_wavelength / 4;
1156 *   if (std::abs(p[0]) >= (layer_center - layer_width / 2) &&
1157 *   std::abs(p[0]) <= (layer_center + layer_width / 2))
1158 *   {
1159 *   return material_b_rho;
1160 *   }
1161 *   }
1162 *  
1163 * @endcode
1164 *
1165 * and finally the default is material_a.
1166 *
1167 * @code
1168 *   return material_a_rho;
1169 *   }
1170 *  
1171 *  
1172 *  
1173 * @endcode
1174 *
1175 *
1176 * <a name="step_62-TheParametersclassimplementation"></a>
1177 * <h4>The `Parameters` class implementation</h4>
1178 *
1179
1180 *
1181 * The constructor reads all the parameters from the HDF5::Group `data` using
1182 * the HDF5::Group::get_attribute() function.
1183 *
1184 * @code
1185 *   template <int dim>
1186 *   Parameters<dim>::Parameters(HDF5::Group &data)
1187 *   : data(data)
1188 *   , simulation_name(data.get_attribute<std::string>("simulation_name"))
1189 *   , save_vtu_files(data.get_attribute<bool>("save_vtu_files"))
1190 *   , start_frequency(data.get_attribute<double>("start_frequency"))
1191 *   , stop_frequency(data.get_attribute<double>("stop_frequency"))
1192 *   , nb_frequency_points(data.get_attribute<int>("nb_frequency_points"))
1193 *   , lambda(data.get_attribute<double>("lambda"))
1194 *   , mu(data.get_attribute<double>("mu"))
1195 *   , dimension_x(data.get_attribute<double>("dimension_x"))
1196 *   , dimension_y(data.get_attribute<double>("dimension_y"))
1197 *   , nb_probe_points(data.get_attribute<int>("nb_probe_points"))
1198 *   , grid_level(data.get_attribute<int>("grid_level"))
1199 *   , probe_start_point(data.get_attribute<double>("probe_pos_x"),
1200 *   data.get_attribute<double>("probe_pos_y") -
1201 *   data.get_attribute<double>("probe_width_y") / 2)
1202 *   , probe_stop_point(data.get_attribute<double>("probe_pos_x"),
1203 *   data.get_attribute<double>("probe_pos_y") +
1204 *   data.get_attribute<double>("probe_width_y") / 2)
1205 *   , right_hand_side(data)
1206 *   , pml(data)
1207 *   , rho(data)
1208 *   {}
1209 *  
1210 *  
1211 *  
1212 * @endcode
1213 *
1214 *
1215 * <a name="step_62-TheQuadratureCacheclassimplementation"></a>
1216 * <h4>The `QuadratureCache` class implementation</h4>
1217 *
1218
1219 *
1220 * We need to reserve enough space for the mass and stiffness matrices and the
1221 * right hand side vector.
1222 *
1223 * @code
1224 *   template <int dim>
1225 *   QuadratureCache<dim>::QuadratureCache(const unsigned int dofs_per_cell)
1226 *   : dofs_per_cell(dofs_per_cell)
1227 *   , mass_coefficient(dofs_per_cell, dofs_per_cell)
1228 *   , stiffness_coefficient(dofs_per_cell, dofs_per_cell)
1229 *   , right_hand_side(dofs_per_cell)
1230 *   {}
1231 *  
1232 *  
1233 *  
1234 * @endcode
1235 *
1236 *
1237 * <a name="step_62-ImplementationoftheElasticWaveclass"></a>
1238 * <h3>Implementation of the `ElasticWave` class</h3>
1239 *
1240
1241 *
1242 *
1243 * <a name="step_62-Constructor"></a>
1244 * <h4>Constructor</h4>
1245 *
1246
1247 *
1248 * This is very similar to the constructor of @ref step_40 "step-40". In addition we create
1249 * the HDF5 datasets `frequency_dataset`, `position_dataset` and
1250 * `displacement`. Note the use of the `template` keyword for the creation of
1251 * the HDF5 datasets. It is a C++ requirement to use the `template` keyword in
1252 * order to treat `create_dataset` as a dependent template name.
1253 *
1254 * @code
1255 *   template <int dim>
1256 *   ElasticWave<dim>::ElasticWave(const Parameters<dim> &parameters)
1257 *   : parameters(parameters)
1258 *   , mpi_communicator(MPI_COMM_WORLD)
1259 *   , triangulation(mpi_communicator,
1260 *   typename Triangulation<dim>::MeshSmoothing(
1261 *   Triangulation<dim>::smoothing_on_refinement |
1262 *   Triangulation<dim>::smoothing_on_coarsening))
1263 *   , quadrature_formula(2)
1264 *   , fe(FE_Q<dim>(1) ^ dim)
1265 *   , dof_handler(triangulation)
1266 *   , frequency(parameters.nb_frequency_points)
1267 *   , probe_positions(parameters.nb_probe_points, dim)
1268 *   , frequency_dataset(parameters.data.template create_dataset<double>(
1269 *   "frequency",
1270 *   std::vector<hsize_t>{parameters.nb_frequency_points}))
1271 *   , probe_positions_dataset(parameters.data.template create_dataset<double>(
1272 *   "position",
1273 *   std::vector<hsize_t>{parameters.nb_probe_points, dim}))
1274 *   , displacement(
1275 *   parameters.data.template create_dataset<std::complex<double>>(
1276 *   "displacement",
1277 *   std::vector<hsize_t>{parameters.nb_probe_points,
1278 *   parameters.nb_frequency_points}))
1279 *   , pcout(std::cout,
1280 *   (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
1281 *   , computing_timer(mpi_communicator,
1282 *   pcout,
1283 *   TimerOutput::never,
1284 *   TimerOutput::wall_times)
1285 *   {}
1286 *  
1287 *  
1288 *  
1289 * @endcode
1290 *
1291 *
1292 * <a name="step_62-ElasticWavesetup_system"></a>
1293 * <h4>ElasticWave::setup_system</h4>
1294 *
1295
1296 *
1297 * There is nothing new in this function, the only difference with @ref step_40 "step-40" is
1298 * that we don't have to apply boundary conditions because we use the PMLs to
1300 *
1301 * @code
1302 *   template <int dim>
1304 *   {
1306 *  
1307 *   dof_handler.distribute_dofs(fe);
1308 *  
1309 *   locally_owned_dofs = dof_handler.locally_owned_dofs();
1312 *  
1313 *   locally_relevant_solution.reinit(locally_owned_dofs,
1315 *   mpi_communicator);
1316 *  
1317 *   system_rhs.reinit(locally_owned_dofs, mpi_communicator);
1318 *  
1319 *   constraints.clear();
1320 *   constraints.reinit(locally_relevant_dofs);
1321 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
1322 *  
1323 *   constraints.close();
1324 *  
1326 *  
1327 *   DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false);
1329 *   locally_owned_dofs,
1330 *   mpi_communicator,
1332 *  
1333 *   system_matrix.reinit(locally_owned_dofs,
1334 *   locally_owned_dofs,
1335 *   dsp,
1336 *   mpi_communicator);
1337 *   }
1338 *  
1339 *  
1340 *  
1341 * @endcode
1342 *
1343 *
1344 * <a name="step_62-ElasticWaveassemble_system"></a>
1346 *
1347
1348 *
1349 * This function is also very similar to @ref step_40 "step-40", though there are notable
1350 * differences. We assemble the system for each frequency/omega step. In the
1351 * first step we set `calculate_quadrature_data = True` and we calculate the
1352 * mass and stiffness matrices and the right hand side vector. In the
1354 *
1355 * @code
1356 *   template <int dim>
1357 *   void ElasticWave<dim>::assemble_system(const double omega,
1358 *   const bool calculate_quadrature_data)
1359 *   {
1360 *   TimerOutput::Scope t(computing_timer, "assembly");
1361 *  
1362 *   FEValues<dim> fe_values(fe,
1366 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1367 *   const unsigned int n_q_points = quadrature_formula.size();
1368 *  
1369 *   FullMatrix<std::complex<double>> cell_matrix(dofs_per_cell, dofs_per_cell);
1370 *   Vector<std::complex<double>> cell_rhs(dofs_per_cell);
1371 *  
1372 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1373 *  
1374 * @endcode
1375 *
1376 * Here we store the value of the right hand side, rho and the PML.
1377 *
1378 * @code
1379 *   std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim));
1380 *   std::vector<double> rho_values(n_q_points);
1381 *   std::vector<Vector<std::complex<double>>> pml_values(
1382 *   n_q_points, Vector<std::complex<double>>(dim));
1383 *  
1384 * @endcode
1385 *
1386 * We calculate the stiffness tensor for the @f$\lambda@f$ and @f$\mu@f$ that have
1389 *
1390 * @code
1392 *   get_stiffness_tensor<dim>(parameters.lambda, parameters.mu);
1393 *  
1394 * @endcode
1395 *
1396 * We use the same method of @ref step_20 "step-20" for vector-valued problems.
1397 *
1398 * @code
1400 *  
1401 *   for (const auto &cell : dof_handler.active_cell_iterators())
1402 *   if (cell->is_locally_owned())
1403 *   {
1404 *   cell_matrix = 0;
1405 *   cell_rhs = 0;
1406 *  
1407 * @endcode
1408 *
1409 * We have to calculate the values of the right hand side, rho and
1410 * the PML only if we are going to calculate the mass and the
1411 * stiffness matrices. Otherwise we can skip this calculation which
1413 *
1414 * @code
1416 *   {
1417 *   fe_values.reinit(cell);
1418 *  
1419 *   parameters.right_hand_side.vector_value_list(
1420 *   fe_values.get_quadrature_points(), rhs_values);
1421 *   parameters.rho.value_list(fe_values.get_quadrature_points(),
1422 *   rho_values);
1423 *   parameters.pml.vector_value_list(
1424 *   fe_values.get_quadrature_points(), pml_values);
1425 *   }
1426 *  
1427 * @endcode
1428 *
1429 * We have done this in @ref step_18 "step-18". Get a pointer to the quadrature
1430 * cache data local to the present cell, and, as a defensive
1431 * measure, make sure that this pointer is within the bounds of the
1432 * global array:
1433 *
1434 * @code
1436 *   reinterpret_cast<QuadratureCache<dim> *>(cell->user_pointer());
1438 *   ExcInternalError());
1440 *   ExcInternalError());
1441 *   for (unsigned int q = 0; q < n_q_points; ++q)
1442 *   {
1443 * @endcode
1444 *
1445 * The quadrature_data variable is used to store the mass and
1446 * stiffness matrices, the right hand side vector and the value
1447 * of `JxW`.
1448 *
1449 * @code
1450 *   QuadratureCache<dim> &quadrature_data =
1452 *  
1453 * @endcode
1454 *
1455 * Below we declare the force vector and the parameters of the
1456 * PML @f$s@f$ and @f$\xi@f$.
1457 *
1458 * @code
1461 *   std::complex<double> xi(1, 0);
1462 *  
1463 * @endcode
1464 *
1466 * step.
1467 *
1468 * @code
1470 *   {
1471 * @endcode
1472 *
1473 * Store the value of `JxW`.
1474 *
1475 * @code
1476 *   quadrature_data.JxW = fe_values.JxW(q);
1477 *  
1478 *   for (unsigned int component = 0; component < dim; ++component)
1479 *   {
1480 * @endcode
1481 *
1482 * Convert vectors to tensors and calculate xi
1483 *
1484 * @code
1485 *   force[component] = rhs_values[q][component];
1486 *   s[component] = pml_values[q][component];
1487 *   xi *= s[component];
1488 *   }
1489 *  
1490 * @endcode
1491 *
1492 * Here we calculate the @f$\alpha_{mnkl}@f$ and @f$\beta_{mnkl}@f$
1493 * tensors.
1494 *
1495 * @code
1498 *   for (unsigned int m = 0; m < dim; ++m)
1499 *   for (unsigned int n = 0; n < dim; ++n)
1500 *   for (unsigned int k = 0; k < dim; ++k)
1501 *   for (unsigned int l = 0; l < dim; ++l)
1502 *   {
1503 *   alpha[m][n][k][l] = xi *
1504 *   stiffness_tensor[m][n][k][l] /
1505 *   (2.0 * s[n] * s[k]);
1506 *   beta[m][n][k][l] = xi *
1507 *   stiffness_tensor[m][n][k][l] /
1508 *   (2.0 * s[n] * s[l]);
1509 *   }
1510 *  
1511 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1512 *   {
1513 *   const Tensor<1, dim> phi_i =
1514 *   fe_values[displacement].value(i, q);
1515 *   const Tensor<2, dim> grad_phi_i =
1516 *   fe_values[displacement].gradient(i, q);
1517 *  
1518 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1519 *   {
1520 *   const Tensor<1, dim> phi_j =
1521 *   fe_values[displacement].value(j, q);
1522 *   const Tensor<2, dim> grad_phi_j =
1523 *   fe_values[displacement].gradient(j, q);
1524 *  
1525 * @endcode
1526 *
1527 * calculate the values of the @ref GlossMassMatrix "mass matrix".
1528 *
1529 * @code
1530 *   quadrature_data.mass_coefficient[i][j] =
1531 *   rho_values[q] * xi * phi_i * phi_j;
1532 *  
1533 * @endcode
1534 *
1535 * Loop over the @f$mnkl@f$ indices of the stiffness
1536 * tensor.
1537 *
1538 * @code
1539 *   std::complex<double> stiffness_coefficient = 0;
1540 *   for (unsigned int m = 0; m < dim; ++m)
1541 *   for (unsigned int n = 0; n < dim; ++n)
1542 *   for (unsigned int k = 0; k < dim; ++k)
1543 *   for (unsigned int l = 0; l < dim; ++l)
1544 *   {
1545 * @endcode
1546 *
1547 * Here we calculate the stiffness matrix.
1550 * gradient function (see the
1551 * [documentation](https://www.dealii.org/current/doxygen/deal.II/group__vector__valued.html))
1552 * which is a <code>Tensor@<2,dim@></code>.
1553 * The matrix @f$G_{ij}@f$ consists of entries
1554 * @f[
1555 * G_{ij}=
1558 * @f]
1559 * Note the position of the indices @f$i@f$ and
1560 * @f$j@f$ and the notation that we use in this
1562 * stiffness tensor is not symmetric, it is
1563 * very easy to make a mistake.
1564 *
1565 * @code
1567 *   grad_phi_i[m][n] *
1568 *   (alpha[m][n][k][l] * grad_phi_j[l][k] +
1569 *   beta[m][n][k][l] * grad_phi_j[k][l]);
1570 *   }
1571 *  
1572 * @endcode
1573 *
1574 * We save the value of the stiffness matrix in
1575 * quadrature_data
1576 *
1577 * @code
1578 *   quadrature_data.stiffness_coefficient[i][j] =
1580 *   }
1581 *  
1582 * @endcode
1583 *
1584 * and the value of the right hand side in
1585 * quadrature_data.
1586 *
1587 * @code
1588 *   quadrature_data.right_hand_side[i] =
1589 *   phi_i * force * fe_values.JxW(q);
1590 *   }
1591 *   }
1592 *  
1593 * @endcode
1594 *
1595 * We loop again over the degrees of freedom of the cells to
1596 * calculate the system matrix. These loops are really quick
1598 * matrices, only the value of @f$\omega@f$ changes.
1599 *
1600 * @code
1601 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1602 *   {
1603 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1604 *   {
1605 *   std::complex<double> matrix_sum = 0;
1606 *   matrix_sum += -Utilities::fixed_power<2>(omega) *
1607 *   quadrature_data.mass_coefficient[i][j];
1608 *   matrix_sum += quadrature_data.stiffness_coefficient[i][j];
1609 *   cell_matrix(i, j) += matrix_sum * quadrature_data.JxW;
1610 *   }
1611 *   cell_rhs(i) += quadrature_data.right_hand_side[i];
1612 *   }
1613 *   }
1614 *   cell->get_dof_indices(local_dof_indices);
1615 *   constraints.distribute_local_to_global(cell_matrix,
1616 *   cell_rhs,
1617 *   local_dof_indices,
1618 *   system_matrix,
1619 *   system_rhs);
1620 *   }
1621 *  
1622 *   system_matrix.compress(VectorOperation::add);
1624 *   }
1625 *  
1626 * @endcode
1627 *
1628 *
1629 * <a name="step_62-ElasticWavesolve"></a>
1630 * <h4>ElasticWave::solve</h4>
1631 *
1632
1633 *
1634 * This is even more simple than in @ref step_40 "step-40". We use the parallel direct solver
1635 * MUMPS which requires less options than an iterative solver. The drawback is
1638 * multigrid method is a well known approach to precondition this system, but
1639 * this is beyond the scope of this tutorial.
1640 *
1641 * @code
1642 *   template <int dim>
1644 *   {
1647 *   locally_owned_dofs, mpi_communicator);
1648 *  
1649 *   SolverControl solver_control;
1650 *   PETScWrappers::SparseDirectMUMPS solver(solver_control, mpi_communicator);
1651 *   solver.solve(system_matrix, completely_distributed_solution, system_rhs);
1652 *  
1653 *   pcout << " Solved in " << solver_control.last_step() << " iterations."
1654 *   << std::endl;
1655 *   constraints.distribute(completely_distributed_solution);
1657 *   }
1658 *  
1659 * @endcode
1660 *
1661 *
1662 * <a name="step_62-ElasticWaveinitialize_position_vector"></a>
1664 *
1665
1666 *
1667 * We use this function to calculate the values of the position vector.
1668 *
1669 * @code
1670 *   template <int dim>
1672 *   {
1673 *   for (unsigned int position_idx = 0;
1674 *   position_idx < parameters.nb_probe_points;
1675 *   ++position_idx)
1676 *   {
1677 * @endcode
1678 *
1679 * Because of the way the operator + and - are overloaded to subtract
1680 * two points, the following has to be done:
1681 * `Point_b<dim> + (-Point_a<dim>)`
1682 *
1683 * @code
1684 *   const Point<dim> p =
1685 *   (position_idx / ((double)(parameters.nb_probe_points - 1))) *
1686 *   (parameters.probe_stop_point + (-parameters.probe_start_point)) +
1687 *   parameters.probe_start_point;
1688 *   probe_positions[position_idx][0] = p[0];
1689 *   probe_positions[position_idx][1] = p[1];
1690 *   if (dim == 3)
1691 *   {
1692 *   probe_positions[position_idx][2] = p[2];
1693 *   }
1694 *   }
1695 *   }
1696 *  
1697 * @endcode
1698 *
1699 *
1700 * <a name="step_62-ElasticWavestore_frequency_step_data"></a>
1702 *
1703
1704 *
1705 * This function stores in the HDF5 file the measured energy by the probe.
1706 *
1707 * @code
1708 *   template <int dim>
1709 *   void
1711 *   {
1712 *   TimerOutput::Scope t(computing_timer, "store_frequency_step_data");
1713 *  
1714 * @endcode
1715 *
1716 * We store the displacement in the @f$x@f$ direction; the displacement in the
1717 * @f$y@f$ direction is negligible.
1718 *
1719 * @code
1720 *   const unsigned int probe_displacement_component = 0;
1721 *  
1722 * @endcode
1723 *
1725 * points of the probe that are located in locally owned cells. The vector
1727 *
1728 * @code
1729 *   std::vector<hsize_t> coordinates;
1730 *   std::vector<std::complex<double>> displacement_data;
1731 *  
1732 *   const auto &mapping = get_default_linear_mapping(triangulation);
1735 *   std::vector<bool> marked_vertices = {};
1736 *   const double tolerance = 1.e-10;
1737 *  
1738 *   for (unsigned int position_idx = 0;
1739 *   position_idx < parameters.nb_probe_points;
1740 *   ++position_idx)
1741 *   {
1743 *   for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1744 *   {
1746 *   }
1747 *   bool point_in_locally_owned_cell = false;
1748 *   {
1750 *   cache, point, cell_hint, marked_vertices, tolerance);
1751 *   if (cell_and_ref_point.first.state() == IteratorState::valid)
1752 *   {
1753 *   cell_hint = cell_and_ref_point.first;
1755 *   cell_and_ref_point.first->is_locally_owned();
1756 *   }
1757 *   }
1759 *   {
1760 * @endcode
1761 *
1762 * Then we can store the values of the displacement in the points of
1764 *
1765 * @code
1767 *   VectorTools::point_value(dof_handler,
1769 *   point,
1770 *   tmp_vector);
1771 *   coordinates.emplace_back(position_idx);
1772 *   coordinates.emplace_back(frequency_idx);
1773 *   displacement_data.emplace_back(
1775 *   }
1776 *   }
1777 *  
1778 * @endcode
1779 *
1780 * We write the displacement data in the HDF5 file. The call
1783 *
1784 * @code
1785 *   if (coordinates.size() > 0)
1786 *   {
1787 *   displacement.write_selection(displacement_data, coordinates);
1788 *   }
1789 * @endcode
1790 *
1791 * Therefore even if the process has no data to write it has to participate
1793 * Note that we have to specify the data type, in this case
1794 * `std::complex<double>`.
1795 *
1796 * @code
1797 *   else
1798 *   {
1799 *   displacement.write_none<std::complex<double>>();
1800 *   }
1801 *  
1802 * @endcode
1803 *
1804 * If the variable `save_vtu_files` in the input file equals `True` then all
1806 * been described in @ref step_40 "step-40".
1807 *
1808 * @code
1809 *   if (parameters.save_vtu_files)
1810 *   {
1811 *   std::vector<std::string> solution_names(dim, "displacement");
1812 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
1815 *  
1816 *   DataOut<dim> data_out;
1817 *   data_out.add_data_vector(dof_handler,
1820 *   interpretation);
1821 *   Vector<float> subdomain(triangulation.n_active_cells());
1822 *   for (unsigned int i = 0; i < subdomain.size(); ++i)
1823 *   subdomain(i) = triangulation.locally_owned_subdomain();
1824 *   data_out.add_data_vector(subdomain, "subdomain");
1825 *  
1826 *   std::vector<Vector<double>> force(
1827 *   dim, Vector<double>(triangulation.n_active_cells()));
1828 *   std::vector<Vector<double>> pml(
1829 *   dim, Vector<double>(triangulation.n_active_cells()));
1830 *   Vector<double> rho(triangulation.n_active_cells());
1831 *  
1832 *   for (auto &cell : triangulation.active_cell_iterators())
1833 *   {
1834 *   if (cell->is_locally_owned())
1835 *   {
1836 *   for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1837 *   {
1838 *   force[dim_idx](cell->active_cell_index()) =
1839 *   parameters.right_hand_side.value(cell->center(), dim_idx);
1840 *   pml[dim_idx](cell->active_cell_index()) =
1841 *   parameters.pml.value(cell->center(), dim_idx).imag();
1842 *   }
1843 *   rho(cell->active_cell_index()) =
1844 *   parameters.rho.value(cell->center());
1845 *   }
1846 * @endcode
1847 *
1848 * And on the cells that we are not interested in, set the
1849 * respective value to a bogus value in order to make sure that if
1851 * looking at the graphical output:
1852 *
1853 * @code
1854 *   else
1855 *   {
1856 *   for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1857 *   {
1858 *   force[dim_idx](cell->active_cell_index()) = -1e+20;
1859 *   pml[dim_idx](cell->active_cell_index()) = -1e+20;
1860 *   }
1861 *   rho(cell->active_cell_index()) = -1e+20;
1862 *   }
1863 *   }
1864 *  
1865 *   for (unsigned int dim_idx = 0; dim_idx < dim; ++dim_idx)
1866 *   {
1867 *   data_out.add_data_vector(force[dim_idx],
1868 *   "force_" + std::to_string(dim_idx));
1869 *   data_out.add_data_vector(pml[dim_idx],
1870 *   "pml_" + std::to_string(dim_idx));
1871 *   }
1872 *   data_out.add_data_vector(rho, "rho");
1873 *  
1874 *   data_out.build_patches();
1875 *  
1876 *   std::stringstream frequency_idx_stream;
1877 *   const unsigned int nb_number_positions =
1878 *   ((unsigned int)std::log10(parameters.nb_frequency_points)) + 1;
1880 *   << std::setfill('0') << frequency_idx;
1881 *   const std::string filename = (parameters.simulation_name + "_" +
1882 *   frequency_idx_stream.str() + ".vtu");
1883 *   data_out.write_vtu_in_parallel(filename, mpi_communicator);
1884 *   }
1885 *   }
1886 *  
1887 *  
1888 *  
1889 * @endcode
1890 *
1891 *
1892 * <a name="step_62-ElasticWaveoutput_results"></a>
1894 *
1895
1896 *
1898 *
1899 * @code
1900 *   template <int dim>
1902 *   {
1903 * @endcode
1904 *
1905 * The vectors `frequency` and `position` are the same for all the
1909 *
1910 * @code
1911 *   if (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)
1912 *   {
1915 *   }
1916 *   else
1917 *   {
1918 *   frequency_dataset.write_none<double>();
1919 *   probe_positions_dataset.write_none<double>();
1920 *   }
1921 *   }
1922 *  
1923 *  
1924 *  
1925 * @endcode
1926 *
1927 *
1928 * <a name="step_62-ElasticWavesetup_quadrature_cache"></a>
1930 *
1931
1932 *
1933 * We use this function at the beginning of our computations to set up initial
1934 * values of the cache variables. This function has been described in @ref step_18 "step-18".
1935 * There are no differences with the function of @ref step_18 "step-18".
1936 *
1937 * @code
1938 *   template <int dim>
1940 *   {
1941 *   triangulation.clear_user_data();
1942 *  
1943 *   {
1944 *   std::vector<QuadratureCache<dim>> tmp;
1945 *   quadrature_cache.swap(tmp);
1946 *   }
1947 *  
1948 *   quadrature_cache.resize(triangulation.n_locally_owned_active_cells() *
1949 *   quadrature_formula.size(),
1950 *   QuadratureCache<dim>(fe.n_dofs_per_cell()));
1951 *   unsigned int cache_index = 0;
1952 *   for (const auto &cell : triangulation.active_cell_iterators())
1953 *   if (cell->is_locally_owned())
1954 *   {
1955 *   cell->set_user_pointer(&quadrature_cache[cache_index]);
1956 *   cache_index += quadrature_formula.size();
1957 *   }
1958 *   Assert(cache_index == quadrature_cache.size(), ExcInternalError());
1959 *   }
1960 *  
1961 *  
1962 *  
1963 * @endcode
1964 *
1965 *
1966 * <a name="step_62-ElasticWavefrequency_sweep"></a>
1968 *
1969
1970 *
1971 * For clarity we divide the function `run` of @ref step_40 "step-40" into the functions
1973 * iteration over the frequency vector.
1974 *
1975 * @code
1976 *   template <int dim>
1978 *   {
1979 *   for (unsigned int frequency_idx = 0;
1980 *   frequency_idx < parameters.nb_frequency_points;
1981 *   ++frequency_idx)
1982 *   {
1983 *   pcout << parameters.simulation_name + " frequency idx: "
1984 *   << frequency_idx << '/' << parameters.nb_frequency_points - 1
1985 *   << std::endl;
1986 *  
1987 *  
1988 *  
1989 *   setup_system();
1990 *   if (frequency_idx == 0)
1991 *   {
1992 *   pcout << " Number of active cells : "
1993 *   << triangulation.n_active_cells() << std::endl;
1994 *   pcout << " Number of degrees of freedom : "
1995 *   << dof_handler.n_dofs() << std::endl;
1996 *   }
1997 *  
1998 *   if (frequency_idx == 0)
1999 *   {
2000 * @endcode
2001 *
2002 * Write the simulation parameters only once
2003 *
2004 * @code
2005 *   parameters.data.set_attribute("active_cells",
2006 *   triangulation.n_active_cells());
2007 *   parameters.data.set_attribute("degrees_of_freedom",
2008 *   dof_handler.n_dofs());
2009 *   }
2010 *  
2011 * @endcode
2012 *
2013 * We calculate the frequency and omega values for this particular step.
2014 *
2015 * @code
2016 *   const double current_loop_frequency =
2017 *   (parameters.start_frequency +
2018 *   frequency_idx *
2019 *   (parameters.stop_frequency - parameters.start_frequency) /
2020 *   (parameters.nb_frequency_points - 1));
2021 *   const double current_loop_omega =
2023 *  
2024 * @endcode
2025 *
2026 * In the first frequency step we calculate the mass and stiffness
2027 * matrices and the right hand side. In the subsequent frequency steps
2029 * time.
2030 *
2031 * @code
2033 *   (frequency_idx == 0) ? true : false);
2034 *   solve();
2035 *  
2038 *  
2039 *   computing_timer.print_summary();
2040 *   computing_timer.reset();
2041 *   pcout << std::endl;
2042 *   }
2043 *   }
2044 *  
2045 *  
2046 *  
2047 * @endcode
2048 *
2049 *
2050 * <a name="step_62-ElasticWaverun"></a>
2051 * <h4>ElasticWave::run</h4>
2052 *
2053
2054 *
2055 * This function is very similar to the one in @ref step_40 "step-40".
2056 *
2057 * @code
2058 *   template <int dim>
2059 *   void ElasticWave<dim>::run()
2060 *   {
2061 *   #ifdef DEBUG
2062 *   pcout << "Debug mode" << std::endl;
2063 *   #else
2064 *   pcout << "Release mode" << std::endl;
2065 *   #endif
2066 *  
2067 *   {
2068 *   Point<dim> p1;
2069 *   p1(0) = -parameters.dimension_x / 2;
2070 *   p1(1) = -parameters.dimension_y / 2;
2071 *   if (dim == 3)
2072 *   {
2073 *   p1(2) = -parameters.dimension_y / 2;
2074 *   }
2075 *   Point<dim> p2;
2076 *   p2(0) = parameters.dimension_x / 2;
2077 *   p2(1) = parameters.dimension_y / 2;
2078 *   if (dim == 3)
2079 *   {
2080 *   p2(2) = parameters.dimension_y / 2;
2081 *   }
2082 *   std::vector<unsigned int> divisions(dim);
2083 *   divisions[0] = int(parameters.dimension_x / parameters.dimension_y);
2084 *   divisions[1] = 1;
2085 *   if (dim == 3)
2086 *   {
2087 *   divisions[2] = 1;
2088 *   }
2090 *   divisions,
2091 *   p1,
2092 *   p2);
2093 *   }
2094 *  
2095 *   triangulation.refine_global(parameters.grid_level);
2096 *  
2098 *  
2100 *  
2101 *   frequency_sweep();
2102 *  
2103 *   output_results();
2104 *   }
2105 *   } // namespace step62
2106 *  
2107 *  
2108 *  
2109 * @endcode
2110 *
2111 *
2112 * <a name="step_62-Themainfunction"></a>
2113 * <h4>The main function</h4>
2114 *
2115
2116 *
2117 * The main function is very similar to the one in @ref step_40 "step-40".
2118 *
2119 * @code
2120 *   int main(int argc, char *argv[])
2121 *   {
2122 *   try
2123 *   {
2124 *   using namespace dealii;
2125 *   const unsigned int dim = 2;
2126 *  
2128 *  
2129 *   HDF5::File data_file("results.h5",
2131 *   MPI_COMM_WORLD);
2132 *   auto data = data_file.create_group("data");
2133 *  
2134 * @endcode
2135 *
2138 *
2139 * @code
2140 *   const std::array<std::string, 2> group_names{
2141 *   {"displacement", "calibration"}};
2142 *   for (const std::string &group_name : group_names)
2143 *   {
2144 * @endcode
2145 *
2146 * For each of these two group names, we now create the group and put
2149 * - The dimensions of the waveguide (in @f$x@f$ and @f$y@f$ directions)
2150 * - The position of the probe (in @f$x@f$ and @f$y@f$ directions)
2151 * - The number of points in the probe
2152 * - The global refinement level
2154 * - The number of mirror pairs
2155 * - The material properties
2156 * - The force parameters
2157 * - The PML parameters
2158 * - The frequency parameters
2159 *
2160
2161 *
2162 *
2163 * @code
2164 *   auto group = data.create_group(group_name);
2165 *  
2166 *   group.set_attribute<double>("dimension_x", 2e-5);
2167 *   group.set_attribute<double>("dimension_y", 2e-8);
2168 *   group.set_attribute<double>("probe_pos_x", 8e-6);
2169 *   group.set_attribute<double>("probe_pos_y", 0);
2170 *   group.set_attribute<double>("probe_width_y", 2e-08);
2171 *   group.set_attribute<unsigned int>("nb_probe_points", 5);
2172 *   group.set_attribute<unsigned int>("grid_level", 1);
2173 *   group.set_attribute<double>("cavity_resonance_frequency", 20e9);
2174 *   group.set_attribute<unsigned int>("nb_mirror_pairs", 15);
2175 *  
2176 *   group.set_attribute<double>("poissons_ratio", 0.27);
2177 *   group.set_attribute<double>("youngs_modulus", 270000000000.0);
2178 *   group.set_attribute<double>("material_a_rho", 3200);
2179 *  
2180 *   if (group_name == "displacement")
2181 *   group.set_attribute<double>("material_b_rho", 2000);
2182 *   else
2183 *   group.set_attribute<double>("material_b_rho", 3200);
2184 *  
2185 *   group.set_attribute(
2186 *   "lambda",
2187 *   group.get_attribute<double>("youngs_modulus") *
2188 *   group.get_attribute<double>("poissons_ratio") /
2189 *   ((1 + group.get_attribute<double>("poissons_ratio")) *
2190 *   (1 - 2 * group.get_attribute<double>("poissons_ratio"))));
2191 *   group.set_attribute("mu",
2192 *   group.get_attribute<double>("youngs_modulus") /
2193 *   (2 * (1 + group.get_attribute<double>(
2194 *   "poissons_ratio"))));
2195 *  
2196 *   group.set_attribute<double>("max_force_amplitude", 1e26);
2197 *   group.set_attribute<double>("force_sigma_x", 1e-7);
2198 *   group.set_attribute<double>("force_sigma_y", 1);
2199 *   group.set_attribute<double>("max_force_width_x", 3e-7);
2200 *   group.set_attribute<double>("max_force_width_y", 2e-8);
2201 *   group.set_attribute<double>("force_x_pos", -8e-6);
2202 *   group.set_attribute<double>("force_y_pos", 0);
2203 *  
2204 *   group.set_attribute<bool>("pml_x", true);
2205 *   group.set_attribute<bool>("pml_y", false);
2206 *   group.set_attribute<double>("pml_width_x", 1.8e-6);
2207 *   group.set_attribute<double>("pml_width_y", 5e-7);
2208 *   group.set_attribute<double>("pml_coeff", 1.6);
2209 *   group.set_attribute<unsigned int>("pml_coeff_degree", 2);
2210 *  
2211 *   group.set_attribute<double>("center_frequency", 20e9);
2212 *   group.set_attribute<double>("frequency_range", 0.5e9);
2213 *   group.set_attribute<double>(
2214 *   "start_frequency",
2215 *   group.get_attribute<double>("center_frequency") -
2216 *   group.get_attribute<double>("frequency_range") / 2);
2217 *   group.set_attribute<double>(
2218 *   "stop_frequency",
2219 *   group.get_attribute<double>("center_frequency") +
2220 *   group.get_attribute<double>("frequency_range") / 2);
2221 *   group.set_attribute<unsigned int>("nb_frequency_points", 400);
2222 *  
2223 *   if (group_name == "displacement")
2224 *   group.set_attribute<std::string>("simulation_name",
2225 *   "phononic_cavity_displacement");
2226 *   else
2227 *   group.set_attribute<std::string>("simulation_name",
2228 *   "phononic_cavity_calibration");
2229 *  
2230 *   group.set_attribute<bool>("save_vtu_files", false);
2231 *   }
2232 *  
2233 *   {
2234 * @endcode
2235 *
2236 * Displacement simulation. The parameters are read from the
2238 * group.
2239 *
2240 * @code
2241 *   auto displacement = data.open_group("displacement");
2243 *  
2245 *   elastic_problem.run();
2246 *   }
2247 *  
2248 *   {
2249 * @endcode
2250 *
2251 * Calibration simulation. The parameters are read from the calibration
2252 * HDF5 group and the results are saved in the same HDF5 group.
2253 *
2254 * @code
2255 *   auto calibration = data.open_group("calibration");
2257 *  
2259 *   elastic_problem.run();
2260 *   }
2261 *   }
2262 *   catch (std::exception &exc)
2263 *   {
2264 *   std::cerr << std::endl
2265 *   << std::endl
2266 *   << "----------------------------------------------------"
2267 *   << std::endl;
2268 *   std::cerr << "Exception on processing: " << std::endl
2269 *   << exc.what() << std::endl
2270 *   << "Aborting!" << std::endl
2271 *   << "----------------------------------------------------"
2272 *   << std::endl;
2273 *  
2274 *   return 1;
2275 *   }
2276 *   catch (...)
2277 *   {
2278 *   std::cerr << std::endl
2279 *   << std::endl
2280 *   << "----------------------------------------------------"
2281 *   << std::endl;
2282 *   std::cerr << "Unknown exception!" << std::endl
2283 *   << "Aborting!" << std::endl
2284 *   << "----------------------------------------------------"
2285 *   << std::endl;
2286 *   return 1;
2287 *   }
2288 *  
2289 *   return 0;
2290 *   }
2291 * @endcode
2292<a name="step_62-Results"></a><h1>Results</h1>
2293
2294
2295<a name="step_62-Resonancefrequencyandbandgap"></a><h3>Resonance frequency and bandgap</h3>
2296
2297
2298The results are analyzed in the
2299[Jupyter Notebook](https://github.com/dealii/dealii/blob/master/examples/step-62/step-62.ipynb)
2301@code{.py}
2302h5_file = h5py.File('results.h5', 'r')
2303data = h5_file['data']
2304
2305# Gaussian function that we use to fit the resonance
2307 omega = 2 * constants.pi * freq
2308 omega_m = 2 * constants.pi * freq_m
2310 return max_amplitude * omega_m**2 * gamma**2 / (((omega_m**2 - omega**2)**2 + gamma**2 * omega**2))
2311
2312frequency = data['displacement']['frequency'][...]
2313# Average the probe points
2314displacement = np.mean(data['displacement']['displacement'], axis=0)
2315calibration_displacement = np.mean(data['calibration']['displacement'], axis=0)
2317reflectivity = (np.abs(np.mean(data['displacement']['displacement'][...]**2, axis=0))/
2318 np.abs(np.mean(data['calibration']['displacement'][...]**2, axis=0)))
2319
2320try:
2328 freq_m = fit_result[0]
2332
2333 fig = plt.figure()
2335 plt.xlabel('frequency (GHz)')
2336 plt.ylabel('amplitude^2 (a.u.)')
2337 plt.title('Transmission\n' + 'freq = ' + "%.7g" % (freq_guess / 1e9) + 'GHz Q = ' + "%.6g" % quality_factor)
2338except:
2339 fig = plt.figure()
2341 plt.xlabel('frequency (GHz)')
2342 plt.ylabel('amplitude^2 (a.u.)')
2343 plt.title('Transmission')
2344
2345fig = plt.figure()
2347plt.xlabel('frequency (GHz)')
2348plt.ylabel('phase (rad)')
2349plt.title('Phase (transmission coefficient)\n')
2350
2351plt.show()
2352h5_file.close()
2353@endcode
2354
2356[resonance frequency](https://en.wikipedia.org/wiki/Resonance) and the
2357[the quality factor](https://en.wikipedia.org/wiki/Q_factor).
2358The quality factor is equal to the ratio between the stored energy in the resonator and the energy
2361[full width at half maximum (FWHM)](https://en.wikipedia.org/wiki/Full_width_at_half_maximum).
2364@f[
2365Q = \frac{f_r}{\Delta f} = \frac{\omega_r}{\Delta \omega} =
23662 \pi \times \frac{\text{energy stored}}{\text{energy dissipated per cycle}}
2367@f]
2368
2370has a gaussian shape
2371@f[
2372a^2 = a_\textrm{max}^2\frac{\omega^2\Gamma^2}{(\omega_r^2-\omega^2)^2+\Gamma^2\omega^2}
2373@f]
2376
2377Given the values we have chosen for the parameters, one could estimate the resonance frequency
2378analytically. Indeed, this is then confirmed by what we get in this program:
2382
2383<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.05.png" height="400" />
2384<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.06.png" height="400" />
2385
2388these sorts of devices for.
2389But it is not quite this easy. In practice, there is really only a "band gap", i.e., the device blocks waves other than
2390the desired one at 20GHz only within a certain frequency range. Indeed, to find out how large this "gap" is within
2392input file. We then obtain the following image:
2393
2394<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.07.png" height="400" />
2395
2399
2400<a name="step_62-Modeprofile"></a><h3>Mode profile</h3>
2401
2402
2408
2409<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.08.png" height="400" />
2410
2411On the other hand, out of resonance all the mechanical energy is
2414at the position @f$x=-8\mu\textrm{m}@f$.
2415
2416<img alt="Phononic superlattice cavity" src="https://www.dealii.org/images/steps/developer/step-62.09.png" height="400" />
2417
2418<a name="step_62-Experimentalapplications"></a><h3>Experimental applications</h3>
2419
2420
2422[quantum optomechanics](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.86.1391).
2424but this code can be used as well to simulate "real world" 3D devices such as
2425[micropillar superlattice cavities](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.060101),
2430
2431
2432<a name="step_62-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2433
2434
2435Instead of setting the parameters in the C++ file we could set the parameters
2436using a Python script and save them in the HDF5 file that we will use for
2437the simulations. Then the deal.II program will read the parameters from the
2438HDF5 file.
2439
2440@code{.py}
2441import numpy as np
2442import h5py
2443import matplotlib.pyplot as plt
2444import subprocess
2445import scipy.constants as constants
2446import scipy.optimize
2447
2448# This considerably reduces the size of the svg data
2449plt.rcParams['svg.fonttype'] = 'none'
2450
2451h5_file = h5py.File('results.h5', 'w')
2452data = h5_file.create_group('data')
2453displacement = data.create_group('displacement')
2454calibration = data.create_group('calibration')
2455
2456# Set the parameters
2457for group in [displacement, calibration]:
2458 # Dimensions of the domain
2459 # The waveguide length is equal to dimension_x
2460 group.attrs['dimension_x'] = 2e-5
2461 # The waveguide width is equal to dimension_y
2462 group.attrs['dimension_y'] = 2e-8
2463
2464 # Position of the probe that we use to measure the flux
2465 group.attrs['probe_pos_x'] = 8e-6
2466 group.attrs['probe_pos_y'] = 0
2467 group.attrs['probe_width_y'] = 2e-08
2468
2469 # Number of points in the probe
2470 group.attrs['nb_probe_points'] = 5
2471
2472 # Global refinement
2473 group.attrs['grid_level'] = 1
2474
2475 # Cavity
2476 group.attrs['cavity_resonance_frequency'] = 20e9
2477 group.attrs['nb_mirror_pairs'] = 15
2478
2479 # Material
2480 group.attrs['poissons_ratio'] = 0.27
2481 group.attrs['youngs_modulus'] = 270000000000.0
2482 group.attrs['material_a_rho'] = 3200
2483 if group == displacement:
2484 group.attrs['material_b_rho'] = 2000
2485 else:
2486 group.attrs['material_b_rho'] = 3200
2487 group.attrs['lambda'] = (group.attrs['youngs_modulus'] * group.attrs['poissons_ratio'] /
2488 ((1 + group.attrs['poissons_ratio']) *
2489 (1 - 2 * group.attrs['poissons_ratio'])))
2490 group.attrs['mu']= (group.attrs['youngs_modulus'] / (2 * (1 + group.attrs['poissons_ratio'])))
2491
2492 # Force
2493 group.attrs['max_force_amplitude'] = 1e26
2494 group.attrs['force_sigma_x'] = 1e-7
2495 group.attrs['force_sigma_y'] = 1
2496 group.attrs['max_force_width_x'] = 3e-7
2497 group.attrs['max_force_width_y'] = 2e-8
2498 group.attrs['force_x_pos'] = -8e-6
2499 group.attrs['force_y_pos'] = 0
2500
2501 # PML
2502 group.attrs['pml_x'] = True
2503 group.attrs['pml_y'] = False
2504 group.attrs['pml_width_x'] = 1.8e-6
2505 group.attrs['pml_width_y'] = 5e-7
2506 group.attrs['pml_coeff'] = 1.6
2507 group.attrs['pml_coeff_degree'] = 2
2508
2509 # Frequency sweep
2510 group.attrs['center_frequency'] = 20e9
2511 group.attrs['frequency_range'] = 0.5e9
2512 group.attrs['start_frequency'] = group.attrs['center_frequency'] - group.attrs['frequency_range'] / 2
2513 group.attrs['stop_frequency'] = group.attrs['center_frequency'] + group.attrs['frequency_range'] / 2
2514 group.attrs['nb_frequency_points'] = 400
2515
2516 # Other parameters
2517 if group == displacement:
2518 group.attrs['simulation_name'] = 'phononic_cavity_displacement'
2519 else:
2520 group.attrs['simulation_name'] = 'phononic_cavity_calibration'
2521 group.attrs['save_vtu_files'] = False
2522
2523h5_file.close()
2524@endcode
2525
2526In order to read the HDF5 parameters we have to use the
2527HDF5::File::FileAccessMode::open flag.
2528@code{.py}
2529 HDF5::File data_file("results.h5",
2532 auto data = data_file.open_group("data");
2533@endcode
2534 *
2535 *
2536<a name="step_62-PlainProg"></a>
2537<h1> The plain program</h1>
2538@include "step-62.cc"
2539*/
void write(const Container &data)
Definition hdf5.h:2003
void write_selection(const Container &data, const std::vector< hsize_t > &coordinates)
Definition hdf5.h:2035
void write_none()
Definition hdf5.h:2193
constexpr void clear()
friend class Tensor
Definition tensor.h:865
Point< 2 > first
Definition grid_out.cc:4629
unsigned int level
Definition grid_out.cc:4632
#define Assert(cond, exc)
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())
Definition loop.h:564
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Definition mapping.cc:316
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
const Event initial
Definition event.cc:68
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
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)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
std::pair< typename MeshType< dim, spacedim >::active_cell_iterator, Point< dim > > find_active_cell_around_point(const Mapping< dim, spacedim > &mapping, const MeshType< dim, spacedim > &mesh, const Point< spacedim > &p, const std::vector< bool > &marked_vertices={}, const double tolerance=1.e-10)
Definition hdf5.h:345
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
constexpr char A
constexpr types::blas_int one
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:74
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Number angle(const Tensor< 1, spacedim, Number > &a, const Tensor< 1, spacedim, Number > &b)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:112
void point_value(const DoFHandler< dim, spacedim > &dof, const VectorType &fe_function, const Point< spacedim, double > &point, Vector< typename VectorType::value_type > &value)
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)
long double gamma(const unsigned int n)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:70
constexpr double PI
Definition numbers.h:241
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation